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Chapter 1 

Introduction 



The state of a physical system is the mathematical description that provides a 
complete information on the system. Its knowledge is equivalent to know the 
result of any possible measurement on the system. In Classical Mechanics it 
is always possible, at least in principle, to devise a procedure made of multi- 
ple measurements which fully recovers the state of the system. In Quantum 
Mechanics, on the contrary, this is not possible, due to the fundamental limita- 
tions related to the Heisenberg uncertainty principle ^ E] and the no-cloning 
theorem 0. In fact, on one hand one cannot perform an arbitrary sequence of 
measurements on a single system without inducing on it a back-action of some 
sort. On the other hand, the no-cloning theorem forbids to create a perfect copy 
of the system without already knowing its state in advance. Thus, there is no 
way out, not even in principle, to infer the quantum state of a single system 
without having some prior knowledge on it 

It is possible to estimate the unknown quantum state of a system when many 
identical copies are available in the same state, so that a different measurement 
can be performed on each copy. A procedure of such kind is called quantum 
tomography. The problem of finding a procedure to determine the state of a 
system from multiple copies was first addressed in 1957 by Fano [S], who called 
quorum a set of observables sufficient for a complete determination of the density 
matrix. However, since for a particle it is difficult to devise concretely measur- 
able observables other than position, momentum and energy, the fundamental 
problem of measuring the quantum state has remained at the level of mere spec- 
ulation up to almost ten years ago, when the issue finally entered the realm of 
experimental physics with the pioneering experiments by Raymer's group [Hj in 
the domain of quantum optics. In quantum optics, in fact, using a balanced 
homodyne detector one has the unique opportunity of measuring all possible 
linear combinations of position and momentum of a harmonic oscillator, which 
here represents a single mode of the electromagnetic field. 

The first technique to reconstruct the density matrix from homodyne mea- 
surements — so called homodyne tomography — originated from the observation 
by Vogel and Risken [7] that the collection of probability distributions achieved 
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by homodyne detection is just the Radon transform of the Wigner function W. 
Therefore, as in classical imaging, by Radon transform inversion one can ob- 
tain W, and then from W the matrix elements of the density operator. This 
first method, however, was affected by uncontrollable approximations, since 
arbitrary smoothing parameters are needed for the inverse Radon transform. 
In Ref. the first exact technique was given for measuring experimentally 
the matrix elements of the density operator in the photon-number representa- 
tion, by simply averaging functions of homodyne data. After that, the method 
was further simplified 0, and the feasibility for non-unit quantum efficiency of 
detectors — above some bounds — was established. 

The exact homodyne method has been implemented experimentally to mea- 
sure the photon statistics of a semiconductor laser and the density matrix 
of a squeezed vacuum The success of optical homodyne tomography has 
then stimulated the development of state-reconstruction procedures for atomic 
beams |12) . the experimental determination of the vibrational state of a molecule 
[T^ , of an ensemble of helium atoms , and of a single ion in a Paul trap ■ 

Through quantum tomography the state is perfectly recovered in the limit 
of infinite number of measurements, while in the practical finite-measurements 
case, one can always estimate the statistical error that affects the reconstruction. 
For infinite dimensions the propagation of statistical errors of the density matrix 
elements make them useless for estimating the ensemble average of unbounded 
operators, and a method for estimating the ensemble average of arbitrary ob- 
servable of the field without using the density matrix elements has been derived 
[T?)] . Further insight on the general method of state reconstruction has lead 
to generalize homodyne tomography to any number of modes |17| . and then to 
extend the tomographic method from the harmonic oscillator to an arbitrary 
quantum system using group theory ^1 El 1201 12] ■ A general data analysis 
method has been designed in order to unbias the estimation procedure from any 
known instrumental noise j20j . Moreover, algorithms have been engineered to 
improve the statistical errors on a given sample of experimental data — the so- 
called adaptive tomography |22 — and then max-likelihood strategies have 
been used that improved dramatically statistical errors, however, at the expense 
of some bias in the infinite dimensional case, and of exponential complexity ver- 
sus for the joint tomography of N quantum systems. The latest technical 
developments derive the general tomographic method from spanning sets of 
operators, the previous group theoretical approaches ^1 El 1201 1^ being just 
a particular case of this general method, where the group representation is just 
a device to find suitable operator "orthogonality" and "completeness" relations 
in the linear algebra of operators. Finally, very recently, a method for tomo- 
graphic estimation of the unknown quantum operation of a quantum device has 
been derived (25) . which uses a single fixed input entangled state, which plays 
the role of all possible input states in quantum parallel on the tested device, 
making finally the method a true "quantum radiography" of the functioning of 
a device. 

In this Review we will give a self-contained and unified derivation of the 
methods of quantum tomography, with examples of applications to different 
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kinds of quantum systems, and with particular focus on quantum optics, where 
also some results from experiments are reexamined. The Review is organized as 
follows. 

In Chapter 2 we introduce the generalized Wigner functions [201 HI] and we 
provide the basic elements of detection theory in quantum optics, by giving the 
description of photodetection, homodyne detection, and heterodyne detection. 
As we will see, heterodyne detection also provides a method for estimating the 
ensemble average of polynomials in the field operators, however, it is unsuitable 
for the density matrix elements in the photon-number representation. The ef- 
fect of non unit quantum efficiency is taken into account for all such detection 
schemes. 

In Chapter 3 we give a brief history of quantum tomography, starting with 
the first proposal of Vogel and Risken 7 as the extension to the domain of quan- 
tum optics of the conventional tomographic imaging. As already mentioned, this 
method indirectly recovers the state of the system through the reconstruction 
of the Wigner function, and is affected by uncontrollable bias. The exact ho- 
modyne tomography method of Ref. [5] (successively simplified in Ref. 0) is 
here presented on the basis of the general tomographic method of spanning sets 
of operators of Ref. As another application of the general method, the 

tomography of spin systems [2H| is provided from the group theoretical method 
of Refs. |18l I19L 1^ . In this chapter we include also further developments to 
improve the method, such as the deconvolution techniques of [20] to correct the 
effects of experimental noise by data processing, and the adaptive tomography 
|22| to reduce the statistical fluctuations of tomographic estimators. 

Chapter 4 is devoted to the evaluation from Ref. ^0] of the expectation 
value of arbitrary operators of a single-mode radiation field via homodyne to- 
mography. Here we also report from Ref. [20] the estimation of the added 
noise with respect to the perfect measurement of field observables, for some 
relevant observables, along with a comparison with the noise that would have 
been obtained using heterodyne detection. 

The generalization of Ref. |17j of homodyne tomography to many modes 
of radiation is reviewed in Chapter 5, where it is shown how tomography of a 
multimode field can be performed by using only a single local oscillator with a 
tunable field mode. Some results of Monte Carlo simulations from Ref. ijLTj are 
also shown for the state that describes light from parametric downconversion. 

Chapter 6 reviews some applications of quantum homodyne tomography to 
perform fundamental test of quantum mechanics. The first is the proposal of 
Ref. [30] to measure the nonclassicality of radiation field. The second is the 
scheme of Ref. [32] to test the state reduction rule using light from parametric 
downconversion. Finally, we review some experimental results about tomogra- 
phy of coherent signals with applications to the estimation of losses introduced 
by simple optical components [62\ . 

Chapter 7 reviews the tomographic method of Ref. J^S] to reconstruct the 
quantum operation of a device, such as an amplifier or a measuring device, using 
a single fixed input entangled state, which plays the role of all possible input 
states in a quantum parallel fashion. 
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Chapter 8 is devoted to the reconstruction technique of Ref. [23 based 
on the maximum hkehhood principle. As mentioned, for infinite dimensions 
this method is necessarily biased, however, it is more suited to the estimation 
of a finite number of parameters, as proposed in Ref. or to the state 

determination in the presence of very low number of experimental data |23) . 
Unfortunately, the algorithm of this method has exponential complexity versus 
the number of quantum systems for a joint tomography of many systems. 

Finally, in Chapter 9 we briefiy review Ref. |34| . showing how quantum to- 
mography could be profitably used as a tool for reconstruction and compression 
in classical imaging. 
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Chapter 2 



Wigner functions and 
elements of detection 
theory 

In this chapter we review some simple formulas from Ref. 35 that connect the 
generalized Wigner functions for s-ordering with the density matrix, and vice- 
versa. These formulas prove very useful for quantum mechanical applications 
as, for example, for connecting master equations with Fokker-Planck equations, 
or for evaluating the quantum state from Monte Carlo simulations of Fokker- 
Planck equations, and finally for studying positivity of the generalized Wigner 
functions in the complex plane. Moreover, as we will show in Chapter 3, the 
first proposal of quantum state reconstruction |7j used the Wigner function as 
an intermediate step. 

In the second part of the chapter, we evaluate the probability distribution of 
the photocurrent of photodetectors, balanced homodyne detectors, and hetero- 
dyne detectors. We show that under suitable limits the respective photocurrents 
provide the measurement of the photon number distribution, of the quadrature, 
and of the complex amplitude of a single mode of the electromagnetic field. 
When the effect of non-unit quantum efficiency is taken into account an addi- 
tional noise affects the measurement, giving a Bernoulli convolution for photo- 
detection, and a Gaussian convolution for homodyne and heterodyne detection. 
Extensive use of the results in this chapter will be made in the next chapters 
devoted to quantum homodyne tomography. 

2.1 Wigner functions 

Since Wigner's pioneering work generalized phase-space techniques have 
proved very useful in various branches of physics '36 . As a method to ex- 
press the density operator in terms of c-number functions, the Wigner functions 
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often lead to considerable simplification of the quantum equations of motion, 
as for example, for transforming master equations in operator form into more 
amenable Fokker-Planck differential equations (see, for example, Ref. [HZl)- 
Using the Wigner function one can express quantum-mechanical expectation 
values in form of averages over the complex plane (the classical phase-space), 
the Wigner function playing the role of a c-number quasi-probability distribu- 
tion, which generally can also have negative values. More precisely, the original 
Wigner function allows to easily evaluate expectations of symmetrically ordered 
products of the field operators, corresponding to the Weyl's quantization pro- 
cedure However, with a slight change of the original definition, one defines 
generalized s-ordered Wigner function Ws{a,a*), as follows ,27j 

W,{a,a*)^ [ ^e"^*-"*^+il^l'Tr[i?(A)p] , (2.1) 

where a* denotes the complex conjugate of a, the integral is performed on 
the complex plane with measure d^A = tiReAdlmA, p represents the density 
operator, and 

D(a) = exp(aa''' - a* a) (2.2) 

denotes the displacement operator, where a and a) ([a, — 1) are the annihila- 
tion and creation operators of the field mode of interest. The Wigner functions 
in Eq. (|2.1|) allow to evaluate s-ordered expectation values of the field operators 
through the following relation 

Tr[:(a^)"a":, \ d^a M^,(a, a*) . (2.3) 

JC 

The particular cases s = —1,0, 1 correspond to anti-normal, symmetrical, and 
normal ordering, respectively. In these cases the generalized Wigner function 
Ws{a,a*) are usually denoted by the following symbols and names 

^Q{a,a*) for s = — 1 "Q function" 

W{a, a*) for s = (usual Wigner function) (2.4) 

P{a,a*) for s = 1 "P function" 

For the normal (s = 1) and anti-normal (s = —1) orderings, the following simple 
relations with the density matrix are well known 

Q{a,a*) = {a\p\a) , (2.5) 

p = / (faP{a,a*) \a){a\ , (2.6) 
JC 

where |a) denotes the customary coherent state \a) = D{a)\0), |0) being the 
vacuum state of the field. Among the three particular representations H2.4() . 
the Q function is positively definite and infinitely differentiable (it actually 
represents the probability distribution for ideal joint measurements of position 
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and momentum of the harmonic oscihator: see Sec. 12 .411 . On the other hand, 
the P function is known to be possibly highly singular, and the only pure states 
for which it is positive are the coherent states Finally, the usual Wigner 
function has the remarkable property of providing the probability distribution 
of the quadratures of the field in form of marginal distribution, namely 

/oo 
dlmaW{ae''<',a*e-''^) = ^(Rea|p|Rea)^ , (2.7) 
-OO 

where \x)<^ denotes the (unnormahzable) eigenstate of the field quadrature 

with real eigenvalue x. Notice that any couple of quadratures X^p, X^^^T^j^ 
is canonically conjugate, namely [^(^,-^(^+,7/2] — V^i and it is equivalent to 
position and momentum of a harmonic oscillator. Usually, negative values of 
the Wigner function are viewed as signature of a non-classical state, the most 
eloquent example being the Schrodinger-cat state |40|. whose Wigner function 
is characterized by rapid oscillations around the origin of the complex plane. 
From Eq. (|2.1(l one can notice that all s-ordered Wigner functions are related 
to each other through Gaussian convolution 

W,{a,a*) = \ SmAfi.n ,f , expf-^^|a-/?p) (2.9) 
Jc ""(s - s) \ s' - s J 

= exp(^^^^^^M^,,(a,a*), {s' > s) . (2.10) 

Equation H2.9|l shows the positivity of the generalized Wigner function for s < 
— 1, as a consequence of the positivity of the Q function. From a qualitative point 
of view, the maximum value of s keeping the generalized Wigner functions as 
positive can be considered as an indication of the classical nature of the physical 
state W. 

An equivalent expression for Ws{a,a*) can be derived as follows Eq. 
H2.1|l can be rewritten as 

Wsia,a*) - TT[pD{a)WsD^ia)] , (2.11) 

where 

f ^e^^W" D{\) . (2.12) 
Through the customary Baker-Campbell-Hausdorff (BCH) formula 

expAexpS = expM + B + ^[A,B]j , (2.13) 
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which holds when [A, [A, B]] — [B, [A, B]] = 0, one writes the displacement in 
normal order, and integrating on arg(A) and |A| one obtains 



'^'-1^*) Vs- 1/ 7r(l-s)Vs-iy 

where we used the normal-ordered forms 

:(a^a)": ^ (at)"a" = a'^aia'^a - 1) . . . (a^a - n + 1) , (2.15) 
and the identity 

.^-xa^a. ^ J2 l_^(at)iai ^ (1 „ 2.)a+a_ (2.16) 

1=0 

The density matrix can be recovered from the generalized Wigner functions 
using the following expression 

p=^ [ d2aW,(a, a*)e-T*7l"l%T^-t f (2.17) 

For the proof of Eq. H2.17|l the reader is referred to Ref. In particular, for 
5 = one has the inverse of the Glauber formula 



p = 2 d^aW{a,a*)D{2a){-)'"', (2.18) 
Jc 

whereas for s = 1 one recovers Eq. (|2.6|l that defines the P function. 



2.2 Photodetection 

Light is revealed by exploiting its interaction with atoms/molecules or electrons 
in a solid, and, essentially, each photon ionizes a single atom or promotes an 
electron to a conduction band, and the resulting charge is then amplified to 
produce a measurable pulse. In practice, however, available photodetectors are 
not ideally counting all photons, and their performances is limited by a non-unit 
quantum efficiency C, namely only a fraction C of the incoming photons lead to 
an electric signal, and ultimately to a count: some photons are either reflected 
from the surface of the detector, or are absorbed without being transformed 
into electric pulses. 

Let us consider a light beam entering a photodetector of quantum efficiency 
^, i.e. a detector that transforms just a fraction ^ of the incoming light pulse 
into electric signal. If the detector is small with respect to the coherence length 
of radiation and its window is open for a time interval T, then the Poissonian 
process of counting gives a probability p{m;T) of revealing m photons that 
writes [H] 



p{m; T) — Tr 



,:KCT^exp[-a(T)T]: 



(2.19) 
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where p is the quantum state of Hght, : : denotes the normal ordering of field 
operators, and I{T) is the beam intensity 

I{T) = e(-) (r, t) ■ e(+) (r, t)dt , (2.20) 

given in terms of the positive (negative) frequency part of the electric field 
operator E^+^(r,t) (E^"'(r,t)). The quantity p{t) = CTr [p/(r)] equals the 
probability of a single count during the time interval {t,t + dt). Let us now 
focus our attention to the case of the radiation field excited in a stationary 
state of a single mode at frequency uj. Eq. H2.19|l can be rewritten as 



p^(m) = Tr 



p : j exp(— 77a' a): 



(2.21) 



where the parameter t] — C,chio/V denotes the overall quantum efficiency of the 
photodetector. Using Eqs. (|2.15|) and (|2.16|l one obtains 

p,{m) = £ Pnn ( ^ ) 77"(1 - VT-"" , (2.22) 

n—m ^ ^ 

where 

P7in = {n\p\n) =p,,^i(n) . (2.23) 

Hence, for unit quantum efficiency a photodetector measures the photon num- 
ber distribution of the state, whereas for non unit quantum efficiency the output 
distribution of counts is given by a Bernoulli convolution of the ideal distribu- 
tion. 

The effects of non unit quantum efficiency on the statistics of a photodetec- 
tor, i.e. Eq. I|2.22(l for the output distribution, can be also described by means 
of a simple model in which the realistic photodetector is replaced with an ideal 
photodetector preceded by a beam splitter of transmissivity t = rj. The reflected 
mode is absorbed, whereas the transmitted mode is photo-detected with unit 
quantum efficiency. In order to obtain the probability of measuring m clicks, 
notice that, apart from trivial phase changes, a beam splitter of transmissivity 
T affects the unitary transformation of fields 



where all field modes are considered at the same frequency. Hence, the output 
mode c hitting the detector is given by the linear combination 



c = - Vl - Tb , (2.25) 
and the probability of counts reads 

Pr{m) ^ Tr [C/rP® |0)(0|C/;t|m)(m| (g) 1] 

= Ep-( ^ )(l-^r""^"- (2-26) 
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Eq. H2.26I) reproduces the probability distribution of Eq. I|2.22|l with t = rj. We 
conclude that a photo-detector of quantum efficiency rj is equivalent to a perfect 
photo-detector preceded by a beam splitter of transmissivity rj which accounts 
for the overall losses of the detection process. 



2.3 Balanced homodyne detection 

The balanced homodyne detector provides the measurement of the quadrature 
of the field X,^ in Eq. (|2.8(l . It was proposed by Yuen and Chan 02], and 
subsequently demonstrated by Abbas, Chan and Yee 03]. 




Figure 2.1: Scheme of the balanced homodyne detector. 

The scheme of a balanced homodyne detector is depicted in Fig. 12.11 The 
signal mode a interferes with a strong laser beam mode 6 in a balanced 50/50 
beam splitter. The mode b is the so-called the local oscillator (LO) mode of the 
detector. It operates at the same frequency of a, and is excited by the laser in a 
strong coherent state \z). Since in all experiments that use homodyne detectors 
the signal and the LO beams are generated by a common source, we assume 
that they have a fixed phase relation. In this case the LO phase provides a 
reference for the quadrature measurement, namely we identify the phase of the 
LO with the phase difference between the two modes. As we will see, by tuning 
(f = aigz we can measure the quadrature X^p at different phases. 

After the beam splitter the two modes are detected by two identical pho- 
todetectors (usually linear avalanche photodiodes), and finally the difference of 
photocurrents at zero frequency is electronically processed and rescaled by 2\z\. 
According to Eqs. H2.24|l . the modes at the output of the 50/50 beam splitter 
(r = 1/2) write 

a — b , a -t- 6 
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hence the difference of photocurrents is given by the following operator 



2\z\ 



2\z\ 



(2.28) 



Let us now proceed to evaluate the probability distribution of the output pho- 
tocurrent / for a generic state p of the signal mode a. In the following treatment 
we will follow Refs. 45, 46 . 

Let us consider the moments generating function of the photocurrent / 

x{\)=Tr[p®\z){z\e'^'] , (2.29) 

which provides the probability distribution of / as the Fourier transform 

P[I) = / -e-'^^x(A) . (2.30) 

Using the BCH formula [iZ||lH| for the SU{2) group, namely 

exp (Ca5t - Ca'h) = e^'''^ (l + \Q^Y^'''^^' ^) e'^''^'^ ^ = ^ tan |^| , (2.31) 

one can write the exponential in Eq. H2.29|l in normal-ordered form with respect 
to mode b as follows 



X(A) 



i tanf TTT-T a 

o \2\z \ J 



X 
2H 



a^a-b''b 



(2.32) 



Since mode 6 is in a coherent state \z) the partial trace over b can be evaluated 
as follows 



X(A) 



A 

2Ui 



i tan( TTT-T I za ' 

a \2\z\ J 





\ /"Ml 








cos — — r 






. \2\z\)_ 





Using now Eq. H2.13|l . one can rewrite Eq. 
to o, namely 



X(A) 



^ / izsin(2^)at 



exp 



A 



-2sin^ ( — ) {a^a + \z 



(2.33) 

in normal order with respect 

^iz' sin(2^)a 



,(2.34) 



In the strong-LO limit z — > oo, only the lowest order terms in X/\z\ are retained, 
a^a is neglected with respect to and Eq. (|2.34|) simplifies as follows 



lim x(A) = ( e*2'="^°* exp 



{exp[iXX^])^ , (2.35) 
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where ^p = argz. The generating function in Eq. 
the POVM 



is then equivalent to 



+00 



2^ 



exp[iA(X^ — x)] = 5{X^p — x) = \x)^^p{x\ , 



(2.36) 



namely the projector on the eigenstate of the quadrature X^p with eigenvalue x. 
In conclusion, the balanced homodyne detector achieves the ideal measurement 
of the quadrature in the strong LO limit. In this limit, the probability 
distribution of the output photocurrent / approaches exactly the probability 
distribution p[x, ip) — ^{x\p\x)^ of the quadrature X^, and this for any state p 
of the signal mode a. 

It is easy to take into account non-unit quantum efficiency at detectors. 
According to Eq. H2.25|l one has the replacements 



c 
d 



u, V vacuum modes 



^d - ^/l~7]V , 

and now the output current is rescaled by 2\z\ri, namely 

+ h.c 



1 

2|d 



2ri 



[u + v) 



(2.37) 
(2.38) 



(2.39) 



where only terms containing the strong LO mode h are retained. The POVM is 
then obtained by replacing 



X^ 



X^ 



2-q 



(liy + v^) 



(2.40) 



in Eq. H2.36|l . with ~ {w^e^"^ + we '^'^)/2, w — u,v, and tracing the vacuum 
modes u and v. One then obtains 



n^(a;) = 



+00 



27r 



+°° dX 



27rA2 



: exp 



27rA2 



+00 



2A2 



dx'e \x')^^{x'\ 



2tt 



(2.41) 



where 



At 



4:7] 



(2.42) 



Thus the POVM, and in turn the probability distribution of the output pho- 
tocurrent, are just the Gaussian convolution of the ideal ones with rms A^ = 
V(l-^)/(4?7). 
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2.4 Heterodyne detection 



Heterodyne detection allows to perform the joint measurement of two conjugated 
quadratures of the field 0^1 ISOI' The scheme of the heterodyne detector is 
depicted in Fig. El 

• cos(cOipt) ► Re Z 

• sin(a)iFt) ► imZ 



0(0)) 

Figure 2.2: Scheme of the heterodyne detector. 

A strong local oscillator at frequency lu in a coherent state \a) hits a beam 
splitter with transmissivity t —^ 1, and with the coherent amplitude a such that 
7 = jal -y/ t(1 — r) is kept constant. If the output photocurrent is sampled at the 
intermediate frequency w/_f, just the field modes a and b at frequency uj ± ujip 
are selected by the detector. Modes a and b are usually referred to as signal 
band and image band modes, respectively. In the strong LO limit, upon tracing 
the LO mode, the output photocurrent I{u}if) rescaled by 7 is equivalent to the 
complex operator 

Z=l^ = a-b\ (2.43) 

7 

where the arbitrary phases of modes have been suitably chosen. The hetero- 
dyne photocurrent Z is a normal operator, equivalent to a couple of commuting 
selfadjoint operators 

Z = RcZ + ilmZ , [Z, Z^\ = [ReZ, ImZ] = . (2.44) 

The POVM of the detector is then given by the orthogonal eigenvectors of Z. 
It is here convenient to introduce the notation of Ref. [3] for vectors in the 
tensor product of Hilbert spaces TL®7i 

\A)) = ^""I") ® I*") = ® = ® ' (2.45) 

nm 

where denotes the transposed operator with respect to some pre-chosen 
orthonormal basis. Eq. (|2.45l) exploits the isomorphism between the Hilbert 
space of the Hilbert-Schmidt operators A,B& HSiTL) with scalar product 



a( (»+(%) 




b(co-coip) 
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{A,B) = Tr[AtB], and the Hilbert space of bipartite vectors \ A)), \B)) en®H, 
where one has {{A\B)) = {A,B). 

Using the above notation it is easy to write the eigenvectors of Z with 
eigenvalue z as -^\D{z))). In fact one has 

oo 

Z\D{z))) = {a-b^){Da{z)^h)\I)) = {Da{z)^h){a^b^ +z)Y,\n)^\n) 
= z{Da{z) ® h)\I)) = z\D{z))) . (2.46) 
The orthogonality of such eigenvectors can be verified through the relation 

{{D{z)\D{z'))) = Tv[D\z)D{z')] = 7tS^^\z - z') , (2.47) 
where 6^'^\a) denotes the Dirac deha function over the complex plane 

exp(7a* - 7*a) . (2.48) 

In conventional heterodyne detection the image band mode is in the vacuum 
state, and one is just interested in measuring the field mode a. In this case we 
can evaluate the POVM upon tracing on mode h. One has 

n(z,z*) = -Tn[\D{z))){{D{z)\Ia®\Q){0\] 

= -D{zm{Q\D\z) = -\z){z\, (2.49) 

TT TT 

namely one obtain the projectors on coherent states. The coherent-state POVM 
provides the optimal joint measurement of conjugated quadratures of the field 
\ b'6\ . In fact, heterodyne detection allows to measure the Q-function in Eq. H2.4|) . 
According to Eq. H2.3|l then it provides the expectation value of anti-normal 
ordered field operator. For a state p the expectation value of any quadrature 
X^p is obtained as 

(X^) = TiipX^] = / — Re(ae-*'^)g(a,a*) . (2.50) 

The price to pay for jointly measuring noncommuting observables is an addi- 
tional noise. The rms fiuctuation is evaluated as follows 

f f^[Re{ae~^n?Qia, «*) " (X^f = (AX^) + 1 , (2.51) 

JC 'f' 4 

where (AX^) is the intrinsic noise, and the additional term is usually referred 
to as "the additional 3dB noise due to the joint measure" [Ml 1551 15^ . 

The effect of non-unit quantum efficiency can be taken into account in anal- 
ogous way as in Sec. 12. 31 for homodyne detection. The heterodyne photocurrent 
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is rescaled by an additional factor rj^^'^, and vacuum modes u and v arc intro- 
duced, thus giving |57] 



Upon tracing over modes u and v, one obtain the POVM 



n^(z,z*) = 



^«(0|.(0|e^(^'i-^*)-^*(^''--)|0)„|0), (2.53) 



TT 



7^(1 - ?y) Jc 7rA2 

The probabiHty distribution is then a Gaussian convolution on the complex 
plane of the ideal probability with rms — {1 — irj)/rj. 

Analogously, the coherent-state POVM for conventional heterodyne detec- 
tion with vacuum image band mode is replaced with 

n.(^,^*) - / ^e"^|.')(/| . (2.54) 

From Eqs. H2.9I) we can equivalently say that the heterodyne detection prob- 
ability density is given by the generalized Wigner function Ws{ot,a*)^ with 
s = 1 — |. Notice that for ?y < 1, the average of functions a^a*™ is related to 
the expectation value of a different ordering of field operators. However, one 
has the relevant identity (23 ES| 



(n,m) / \ / \ /j- _ \k 

ia^a-^-.s = E fc! P (^^j ■.iar-'a-^-'-.t , (2.55) 



fc=0 

where (n,m) = min(n, m), and then 



d^aW^i_2(a,a*)a"a*" 

V 



Notice that the measure of the Q-function (or any smoothed version for < 1) 
does not allow to recover the expectation value of any operator through an 
average over heterodyne outcomes. In fact, one needs the admissibility of anti- 
normal ordered expansion and the convergence of the integral in Eq. (|2.56|) . 
In particular, the matrix elements of the density operator cannot be recovered. 
For some operators in which heterodyne measurement is allowed, a comparison 
with quantum homodyne tomography will be given in Sec. 14.31 
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Finally, it is worth mentioning that all results of this section are valid also 
for an image-band mode with the same frequency of the signal. In this case a 
measurement scheme based on multiport homodyne detection should be used 

[SDlEHlEniinilESESllMIESlEni- 
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Chapter 3 

General tomographic 
method 

In this chapter we review the general tomographic method of spanning sets of 
operators of Ref. 24 , and re-derive in this general framework the exact homo- 
dyne tomography method of Ref. "S". In the first section, we first give a brief 
history of quantum tomography, starting with the original proposal of Vogel and 
Risken |Zj , that extended the conventional tomographic imaging to the domain 
of quantum optics. Here we will briefly sketch the conventional imaging tomo- 
graphic method, and show the analogy with the method of Ref. [J- The latter 
achieves the quantum state via the Wigner function, which in turn is obtained 
by inverse Radon transform of the homodyne probability distributions for vary- 
ing phase with respect to the LO. As already mentioned, the Radon transform 
inversion is affected by uncontrollable bias: such limitations and the intrinsic 
unreliability of this method are thoroughly explained in the same section. 

As opposite to the Radon transform method, the first exact method of Ref. 
|S] (and successively refined in Ref. [H]) allows the reconstruction of the density 
matrix p, bypassing the step of the Wigner function, and achieving the matrix 
elements of p — or the expectation of any arbitrary operator — by just averaging 
the pertaining estimators (also called Kernel functions or pattern functions), 
evaluated on the experimental homodyne data. This method will be re-derived 
in Subsec. 13.3.31 as a special case of the general tomographic method of Ref. 23] 
here reviewed in Sect. 13.31 where we introduce the concept of "quorum" , which 
is the complete set of observables whose measurement provides the expectation 
value of any desired operator. Here we also show how some "orthogonality" 
and "completeness" relations in the linear algebra of operators are sufficient 
to individuate a quorum. As another application of the general method, in 
Subsect. 13.3.51 the tomography of spin systems |2H| is reviewed, which was 
originally derived from the group theoretical methods of Refs. ^1 El I2()j . 
Another application is the quantum tomography of a free particle state, given 
in Subsect. 13.3.61 
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In this same chapter, in Sec. 13.41 we include some further developments to 
improve the tomographic method, such as the deconvolution techniques of Ref. 
| 2U| to correct the imperfections of detectors and experimental apparatus with 
a suitable data processing, and the adaptive tomography of Ref. [221 to reduce 
the statistical fluctuations of tomographic estimators, by adapting the averaged 
estimators to the given sample of experimental data. 

The other relevant topics of homodyning observables, multimode tomogra- 
phy, and tomography of quantum operations will be given a separate treatment 
in the following chapters of the Review. 

3.1 Brief historical excursus 

The problem of quantum state determination through repeated measurements 
on identically prepared systems was originally stated by Fano in 1957 jS], who 
first recognized the need of measuring more that two noncommuting observables 
to achieve such purpose. However, it was only with the proposal by Vogel and 
Risken |7] that quantum tomography was born. The first experiments, which 
already showed reconstructions of coherent and squeezed states were performed 
by Michael Raymer's and his group at the University of Oregon ^ . The main 
idea at the basis of the first proposal is that it is possible to extend to the quan- 
tum domain the algorithms that are conventionally used in medical imaging to 
recover two dimensional (mass) distributions from unidimensional projections 
in different directions. This first tomographic method, however, was unreliable 
for the reconstruction of an unknown quantum state, since arbitrary smooth- 
ing parameters were needed in the Radon-transform based imaging procedure. 
The first exact unbiased tomographic method was proposed in Ref. and 
successively simplified in Ref. [Hj. Since then, the new exact method has been 
practically implemented in many experiments, such as the measurement of the 
photon statistics of a semiconductor laser JOI) and the reconstruction of the 
density matrix of a squeezed vacuum ^Jj. The success of optical homodyne 
tomography has then stimulated the development of state-reconstruction proce- 
dures in other quantum harmonic oscillator systems, such as for atomic beams 
|12|. and the vibrational state of a molecule of an ensemble of helium atoms 
| 14| . and of a single ion in a Paul trap |15| . 

After the original exact method, quantum tomography has been general- 
ized to the estimation of arbitrary observable of the field to any num- 
ber of modes ^7], and, finally, to arbitrary quantum systems via group theory 
I19L 1^ I21| . with further improvements such as noise deconvolution |2U| . 
adaptive tomographic methods |23, and the use of max-likelihood strategies 
|2H| . which has made possible to reduce dramatically the number of experimen- 
tal data, up to a factor 10^-^ 10^, with negligible bias for most practical cases of 
interest. Finally, more recently, a method for tomographic estimation of the un- 
known quantum operation of a quantum device has been proposed |25| . where a 
fixed input entangled state plays the role of all input states in a sort of quantum 
parallel fashion. Moreover, as another manifestation of such a quantum paral- 
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lelism, one can also estimate the ensemble average of all operators by measuring 
only one fixed "universal" observable on an extended Hilbert space in a sort 
of quantum hologram |67| . This latest development is based on the general to- 
mographic method of Ref. |24) , where the tomographic reconstruction is based 
on the existence of spanning sets of operators, of which the irreducible unitary 
group representations of the group-methods of Refs. ^1 El HOI 1^ are just a 
special case. 



3.2 Conventional tomographic imaging 

In conventional medical tomography, one collects data in the form of marginal 
distributions of the mass function m{x,y). In the complex plane the marginal 
r(a;, (p) is a projection of the complex function y) on the direction indicated 
by the angle ip G [0,7r], namely 

r{x,ip) — / — m (^{x + iy)e^'^ , {x — iy)e^^'^) . (3.1) 

J — oo ^ 

The collection of marginals for different cp is called "Radon transform". The 
tomography process essentially consists in the inversion of the Radon transform 
H^.lfl . in order to recover the mass function m,{x, y) from the marginals r{x, ip). 
Here we derive inversion of Eq. (|3.1|l . Consider the identity 

m(a,a*)= / (f (3 S'-^^a ~ f3) m{f3, (3*) , (3.2) 



where S^'^\a) denotes the Dirac delta function of Eq. H2.48|l . and m(a, a*) = 
m{x, y), with a = x + iy and a* — x — iy. It is convenient to rewrite Eq. H2.48|l 
as follows 

^(^)(o) ^ r r ^ ^-'^^ - r ? i^i r ^ e---, (3.3) 







4 ./n TT^ 4 



with = Re(a e "^) = — ai^+jr. Then, from Eqs. H3.2|l and l|3.3|l the inverse 
Radon transform is obtained as follows 

m{x,y)^ ^ dx'r{x',^) e^'^C-'"".) . (3.4) 



/ TT I \ 1 r / I , 

'0 J —oo J —oo ^ 

Eq. H3.4|l is conventionally written as 



r dip f^°° 

m(x,y) ^ I — / dx' r{x' ,p) K{x' - a^), (3.5) 

Jo J-oo 

where K{x) is given by 

K{x)= ^\k\e'''^ = -Re dkke'''^ ^--V— , (3.6) 

J—^, 4 2 ./n 2 x 
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with V denoting the Cauchy principal value. Integrating by parts Eq. H3.5II one 
obtains the tomographic formula that is usually found in medical imaging, i.e. 

1 r f^°° 1 d 

m{x,y) = T^- / d-^ P / dx' — — r{x' ,if) , (3.7) 

27r Jq x' - dx' 

which allows the reconstruction of the mass distribution m{x,y) from its pro- 
jections along different directions r(x,ip). 

3.2.1 Extension to the quantum domain 

In the "quantum imaging" process the goal is to reconstruct a quantum state in 
the form of its Wigner function starting from its marginal probability distribu- 
tions. As shown in Sec. 12.11 the Wigner function is a real normalized function 
that is in one-to-one correspondence with the state density operator p. As 
noticed in Eq. 1)2. 7|l . the probability distributions of the quadrature operators 
Xip = (a^e*'^-|-ae~"^)/2 are the marginal probabilities of the Wigner function for 
the state p. Thus, by applying the same procedure outlined in the previous sub- 
section, Vogel and Risken |7j proposed a method to recover the Wigner function 
via an inverse Radon transform from the quadrature probability distributions 
p{x, if), namely 

din /■+°° /--l-oo 17 

W{x,y)^ dx'p{x',ip) _ g^fe(^-xcosv-ysinv) _ (gg) 

Jo ^ J-oo J-oo 4 

[Surprisingly, in the original paper the connection to the tomographic imag- 
ing method was never mentioned]. As shown in Sec. 12.31 the experimental 
measurement of the quadratures of the field is obtained using the homodyne 
detector. The method proposed by Vogel and Risken, namely the inversion of 
the Radon transform, was the one which has been used in the first experiments 

m 

This first method is, however, not reliable for the reconstruction of an un- 
known quantum state, due to the intrinsic unavoidable systematic error related 
to the fact that the integral on k in Eq. 1)3. 8|l is unbounded. In fact, in order 
to evaluate the inverse Radon transform, one would need the analytical form 
of the marginal distribution of the quadrature p{x, if), which, in turn, can only 
be obtained by collecting the experimental data into histograms, and thence 
"spline-ing" them. This, of course, is not an unbiased procedure since the de- 
gree of spline-ing, the width and the number of the histogram bins, and finally 
the number of different phases used to collect the experimental data sample in- 
troduce systematic errors if they are not set above some minimal values, which 
actually depend on the unknown quantum state that one wants to reconstruct. 
Typically, an over-spline-ing will wash-out the quantum features of the state, 
whereas, vice-versa, an under-spline-ing will create negative photon probabilities 
in the reconstruction (see Ref. |8j for details). 

A new exact method was then proposed in Ref. alternative to the Radon- 
transform technique. This approach, referred to as quantum homodyne tomog- 
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raphy, allows to recover the quantum state of the field p — along with any en- 
semble average of arbitrary operators — by directly averaging functions of the 
homodyne data, abolishing the intermediate step of the Wigner function, which 
is the source of all systematic errors. Only statistical errors are present, and 
they can be reduced arbitrarily by collecting more experimental data. This 
exact method will be re-derived from the general tomographic theory in Sec. 



3.3 General method of quantum tomography 

In this section the general method of quantum tomography is explained in de- 
tail. First, we give the basics of Monte Carlo integral theory which are needed 
to implement the tomographic algorithms in actual experiments and in numer- 
ical simulations. Then, we derive the formulas on which all schemes of state 
reconstruction are based. 

3.3.1 Basic statistics 

The aim of quantum tomography is to estimate, for arbitrary quantum system, 
the mean value (O) of a system operator O using only the results of the measure- 
ments on a set of observables {Qx, A S A}, called "quorum". The procedure by 
which this can be obtained needs the estimator or "Kernel function" TZ[0]{x, A) 
which is a function of the eigenvalues x of the quorum operators. Integrating 
the estimator with the probability p(x, A) of having outcome x when measuring 
Qx, the mean value of O is obtained as follows 



where the first integral is performed on the values of A that designate all quorum 
observables, and the second on the eigenvalues of the quorum observable Qx 
determined by the A variable of the outer integral. For discrete set A and/or 
discrete spectrum of the quorum, both integrals in H3.9|l can suitably replaced 
by sums. 

The algorithm to estimate (O) with Eq. (|3.9|l is the following. One chooses 
a quorum operator Qx by drawing A with uniform probability in A and performs 
a measurement, obtaining the result Xi. By repeating the procedure N times, 
one collects the set of experimental data {{xi,Xi), with i = 1, ■ ■ ■ ^ N}, where 
Aj identifies the quorum observable used for the ith measurement, and Xi its 
result. From the same set of data the mean value of any operator O can be 
obtained. In fact, one evaluates the estimator of (O) and the quorum Qx, and 
then samples the double integral of H3.9|) using the limit 



EH 



(O) ^ f dX f dfixix) p{x, A) n[0]ix, A) , 



(3.9) 




(3.10) 
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Of course the finite sum 

1 ^ 

FN^J;^Y.'^^0]{X^,\^)■ (3.11) 

1=1 

gives an approximation of (O) . To estimate the error in the approximation one 
apphes the central hmit theorem that we recaU here. 

Central limit theorem. Consider N statisticahy uncorrelated random 
variables {zi, i = 1,---,A^}, with mean values n^Zi), variances (J^(zi) and 
bounded third order moments. If the variances cr^(zi) are all of the same order 
then the statistical variable "average" y defined as 

1 ^ 



VN 

has mean and variance 



N ■ 

i=l 



1 ^ 1 ^ 

1=1 1=1 

The distribution of j/at approaches asymptotically a Gaussian for iV — > oo. In 
practical cases, the distribution of y can be considered Gaussian already for N 
as low as ~ 10. 

For our needs the hypotheses are met if the estimator TZ[0](xi, Xi) in Eq. 
H3.11|l has limited moments up to the third order, since, even though Xi have 
different probability densities depending on A^, nevertheless, since Xi is also 
random all Zi here given by 

z, = n[0]{x,,X,) (3.14) 
fi{z^) = (O) (3.15) 



have common mean 



and variance 

cr2(z,)= / dX I d^ix{x) pix,X)n''[0]{x,X) - {Oy^. (3.16) 



A 



Using the central limit theorem, we can conclude that the experimental average 
y = Fn in Eq. (|3.11() is a statistical variable distributed as a Gaussian with 
mean value /j.(?/jv) = fJ'{zi) and variance a'^{yN) = ■^cr^(zi). Then the tomo- 
graphic estimation converges with statistical error that decreases as N^^^^. A 
statistically precise estimate of the confidence interval is given by 



'^^V N{N-1) ' (^-^^^ 

with Zi given by Eq. H3.14|l and yN by Eq. (|3.12|l . In order to test that 
the confidence intervals are estimated correctly, one can check that the Fn 
distribution is actually Gaussian. This can be done by comparing the histogram 
of the block data with a Gaussian, or by using the test. 
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3.3.2 Characterization of the quorum 

Different estimation techniques have been proposed tailored to different quan- 
tum systems, such as the radiation field [SI EI) trapped ions and molecular 
vibrational states j68|. spin systems |69| . etc. All the known quantum estima- 
tion techniques can be embodied in the following approach. 

The tomographic reconstruction of an operator O is possible when there 
exists a resolution of the form 

O ^ f d\Tv[OB\\)\C{\) , (3.18) 

where A is a (possibly multidimensional) parameter living on a (discrete or 
continuous) manifold A. The only hypothesis in (|3.18|) is the existence of the 
trace. If, for example, O is a trace-class operator, then we do not need to 
require B{\) to be of Hilbert-Schmidt class, since it is sufficient to require 
B{X) bounded. The operators C(A) are functions of the quorum of observables 
measured for the reconstruction, whereas the operators B{\) form the dual basis 
of the set C(A). The term 

£[0]{\) = Tr [OB^X)] C(A) (3.19) 

represents the quantum estimator for the operator O. The expectation value of 
O is given by the ensemble average 

(O) EE Tr [Op] = / dA Tr [OB^X]] Tr [C(A)p] = [ dX [0](A)) , (3.20) 
Ja J a 

where p is the density matrix of the quantum system under investigation. Notice 
that the quantity Tr [C(A)/9] depends only on the quantum state, and it is related 
to the probability distribution of the measurement outcomes, whereas the term 
Tr [OB'''(A)] depends only on the quantity to be measured. In particular, the 
tomography of the quantum state of a system corresponds to writing Eq. H3.18|l 
for the operators O = \k){n\, {|n)} being a given Hilbert space basis. For a 
given system, the existence of a set of operators C(A), together with its dual 
basis B(X) allows universal quantum estimation, i. e. the reconstruction of any 
operator. 

We now give two characterizations of the sets B{X) and C(A) that are nec- 
essary and sufficient conditions for writing Eq. H3.18|l . 

Condition 1: bi-orthogonality 
Let us consider a complete orthonormal basis of vectors \n) [n = 0, 1,---). 
Formula H3.18|l is equivalent to the bi-orthogonality condition 

/ dX {q\B^X)\p) {m\C{X)\l) = S„,p6ig , (3.21) 
Ja 

where 6ij is the Kronecker delta. Eq. H3.21|l can be straightforwardly generalized 
to a continuous basis. 

Condition 2: completeness 
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If the set of operators C(A) is complete, namely if any operator can be written 
as a linear combination of the C(A) as 



0= dX a(A) C(A) , (3.22) 
Ja 

then Eq. ()3.18fl is also equivalent to the trace condition 

TT[B^X)C{fi)] = S{X,^l) , (3.23) 

where i5(A, /x) is a reproducing kernel for the set -B(A), namely it is a function or 
a tempered distribution which satisfies 



d\ B{X) 6{X, ^) = B{fj,) . (3.24) 

A 

An analogous identity holds for the set of C(A) 

dXC{X) S{X,fi) ^C{n) . (3.25) 

A 

The proofs are straightforward. The completeness condition on the operators 
C(A) is essential for the equivalence of (|3.18|l and H3.23|l . A simple counterex- 
ample is provided by the set of projectors P(A) = |A)(A| over the eigenstates of a 
self-adjoint operator L. In fact, Eq. (|3.23|) is satisfied by C(A) = B{X) = P{X). 
However, since they do not form a complete set in the sense of Eq. 13.22|l . it is 
not possible to express a generic operator in the form X — Jj^dX (A|0|A) |A)(A|. 
If either the set B{X) or the set C(A) satisfy the additional trace condition 

Tr[Bt(^)B(A)] =5(A,/i) , (3.26) 
Tr[cHp)C{X)]=SiX,^l) , (3.27) 

then we have C(A) = B{X) (notice that neither B{X) nor C(A) need to be 
unitary). In this case, Eq. 13.18|l can be rewritten as 

0= [ dXTv [OC^(A)] C(A) . (3.28) 

J A 

A certain number of observables Q\ constitute a quorum when there are 
functions f\{Qx) = C'(A) such that C(A) form an irreducible set. The quantum 
estimator for O in Eq. (|3.19() then writes as a function of the quorum operators 

£[0]{X) = £x[0]{Qx) . (3.29) 

Notice that if a set of observables Q\ constitutes a quorum, than the set of pro- 
jectors \q)\x{q\ over their eigenvectors provides a quorum too, with the measure 
dX in Eq. H3.18I) including the measure d^\{q). Notice also that, even once the 
quorum has been fixed, the unbiased estimator for an operator O will not in 
general be unique, since there can exist functions M{Q\) that satisfies 

dXMiQx) = , (3.30) 

A 
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and that will be called 'null estimators'. Two unbiased estimators that differ 
by a null estimator yield the same results when estimating the operator mean 
value. We will see in Sec. 13.4.21 how the null estimators can be used to reduce 
the statistical noise. 

In terms of the quorum observables Qx Eq. H3.20|l rewrites 

(O) = / dATr[OBt(A)]Tr[A(QA)p] 

= J^dxJdfix{q)p{q,X)Tr[OB\X)]h{q), (3.31) 

where p{q,X) — x{q\p\q)\ is the probability density of getting the outcome q 
from the measurement of Q\ on the state p. Eq. (|3.31|) is equivalent to the 
expression (|3.9|l . with estimator 

7^[0]((?,A) =Tr[OSt(A)] h{q) . (3.32) 

Of course it is of interest to connect a quorum of observables to a resolution 
of the form H3.18|l . since only in this case there can be a feasible reconstruc- 
tion scheme. If a resolution formula is written in terms of a set of selfadjoint 
operators, the set itself constitutes the desired quorum. However, in general a 
quorum of observables is functionally connected to the corresponding resolution 
formula. If the operators C (A) are unitary, then they can always be taken as the 
exponential map of a set of selfadjoint operators, which then are identified with 
our quorum Qx- The quantity Tr[C(A)/c<] is thus connected with the moment 
generating function of the set Qx, and hence to the probability density p{q, A) of 
the measurement outcomes, which play the role of the Radon transform in the 
quantum tomography of the harmonic oscillator. In general, the operators C(A) 
can be any function (neither self-adjoint nor unitary) of observables and, even 
more generally, they may be connected to POVMs rather than observables. 

The dual set B{X) can be obtained from the set C(A) by solving Eq. (|3.23|l . 
For finite quorums, this resorts to a matrix inversion. An alternative procedure 
uses the Gram-Schmidt orthogonalization procedure No such a general 

procedure exists for a continuous spanning set. Many cases, however, satisfy 
conditions (|3.26() and H3.27|l . and thus we can write B{X) — C(A)^. 

3.3.3 Quantum estimation for harmonic-oscillator systems 

The harmonic oscillator models several systems of interest in quantum mechan- 
ics, as the vibrational states of molecules, the motion of an ion in a Paul trap, 
and a single mode radiation field. Different proposals have been suggested in 
order to reconstruct the quantum state of a harmonic system, which all fit the 
framework of the previous section, which is also useful for devising novel estima- 
tion techniques. Here, the basic resolution formula involves the set of displace- 
ment operators D{a) = exp(Q!a^ — a*a), which can be viewed as exponentials of 
the field-quadrature operators — [a) e^'^ + ae~"^)/2. We have shown in Sec. 
12.31 that for a single-mode radiation field X,^ is measured through homodyne 
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detection. For the vibrational tomography of a molecule or a trapped ion X^p 
corresponds to a time-evolved position or momentum. The set of displacement 
operators satisfies Eqs. (|3.23ll and (|3.27|) . since 

T:i[D{a)D\l3)] = 7r^(2)(Q, _ ^) ^ (3 33) 

whereas Eq. H3.28|l reduces to the Glauber formula 

0= / — Tr \OD\a)\ D{a) . (3.34) 

Changing to polar variables a — (— i/2)fce*'^, Eq. H3.34() becomes 

= /•+°° dim ^^.JQ ^^kX,^ ^-^kX, ^ (3 35) 



which shows explicitly the dependence on the quorum X^p. Taking the ensemble 
average of both members and evaluating the trace over the set of eigenvectors 



of X^ , one obtains 



(0)=/ — / dxp{x,^)TZ[0]{x,ip) , (3.36) 



where p{x;ip) — ^{x\p\x)p is the probability distribution of quadratures out- 
come. The estimator of the operator ensemble average (O) is given by 

n[0] {x, ^) = Tr[OK{Xp ~ x)] , (3.37) 

where K{x) is the same as in Eq. (|3.6|l . 

Eq. H3.36|l is the basis of quantum homodyne tomography. Notice that even 
though K{x) is unbounded, however, the matrix element {4'\K{X^ — x)\(f>) can 
be bounded, whence it can be used to sample the matrix element (iplplcj)) of 
the state p, which, according to Sec. 13.3.11 is directly obtained by averaging 
the estimator (I3.37|) over homodyne experimental values. In fact, for bounded 
{il;\K{Xp — x)\<j>), the central limit theorem guarantees that 

{i:\p\cj)) = / ^/ dxp{x,^) {ij\K{Xp-xM (3.38) 

Jo I" 

1 N 

]^ E(V^I^(^^" - , (3.39) 



71 = 



where a;„ is the homodyne outcome measured at phase (p„ and distributed with 
probability p{x,ip). Systematic errors are eliminated by choosing randomly 
each phase ipn at which homodyne measurement is performed. As shown in 
Sec. 13.3.11 for finite number of measurements N, the estimate (|3.39|) of the 
integral in Eq. H3.38|l is Gaussian distributed around the true value 
with statistical error decreasing as N~^^'^. Notice that the measurability of the 
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density operator matrix element depends only on the boundedness of the matrix 
element of the estimator, and that no adjustable parameters are needed in the 
procedure, which thus is unbiased. 

The general procedure for noise deconvolution is presented in Sec. 13.4.11 
However, we give here the main result for the density matrix reconstruction. 
As shown in Sec. 12.31 the effect of the efficiency in homodyne detectors is a 
Gaussian convolution of the ideal probability p{x, ip), as 



p^{x, p) = y ^(^^ dx' e~^^^~^'^\{x', . (3.40) 

The tomographic reconstruction procedure still holds upon replacing p{x^ Lp) 
with Pr,{x, if), so that 



dip 

where now the estimator is 



r dip f^°° 

P = I — / dx pjj{x,ip) Kr^{X^ - x), (3.41) 

Jo ^ J-ao 



K„{x) = -Rc / k dk 6^'='+''=^ . (3.42) 
2 Jo 

In fact, by taking the Fourier transform of both members of Eq. (|3.4l)|) . one can 
easily check that 

d^ '■+°" 



da; Pn{x, Lp) KJX^p - x) 

' 

+ 00 

dxp[x,p) K{X^-x) . (3.43) 



dip 







TT 



Notice that the anti-Gaussian in Eq. (|3.42() causes a much slower convergence 
of the Monte Carlo integral (|3.41l) : the statistical fluctuation will increase ex- 
ponentially for decreasing detector efficiency 77. In order to achieve good recon- 
structions with non- ideal detectors, then one has to collect a larger number of 
data. 

It is clear from Eq. (|3.39|) that the measurability of the density matrix 
depends on the chosen representation and on the quantum efficiency of the 
detectors. For example, for the reconstruction of the density matrix in the Fock 
basis the estimators are given by 

■Rn[\n){n + d\]{x,v)= /"^°° e^fc'-'fc-(„ + rfje^^^^ |n) (3.44) 

J —OC ^ 

/ I r-\-oo 

= e^di^+^) / / dk \k\e'^''"-'"'^k^Li{k') , 
V [n + dy. j_oo 

where L^{x) denotes the generalized Laguerre polynomials. Notice that the 
estimator is bounded only for 77 > 1/2, and below the method would give un- 
bounded statistical errors. However, this bound is well below the values that are 
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reasonably achieved in the lab, where actual homodyne detectors have efficiency 
ranging between 70% and 90% |11[I70| . Moreover, a more efficient algorithm is 
available, that uses the factorization formulas that hold for ry = 1 |^ [72] 

n[\n)(n + d\]{x,^) = e'^^[4xu^{x)vn+dix) (3.45) 
-2^/n + lun+i{x)v„+d{x) - 2Vn+irfTun{x)vn+d+i{x)] , 

where Uj (x) and Vj (x) are the normalizable and unnormalizable eigenfunctions 
of the harmonic oscillator with eigenvalue j, respectively. The noise from quan- 
tum efficiency can be unbiased via the inversion of the Bernoulli convolution, 
which holds for n > 1/2 [THj. 

The use of Eq. (|3.36|l to estimate arbitrary operators through homodyne 
tomography will be the subject of the following Chapter. Notice that Eq. H3.34|l 
cannot be used for unbounded operators, however the estimators also for some 
unbounded operators will be derived in Sec. 14.11 



3.3.4 Some generalizations 

Using condition (|3.23|) one can sec that the Glauber formula can be generalized 
to 



0=1 — TT[0F,Dia)F2]F-'DHa)F-\ 
Jc 



(3.46) 



where Fi and F2 are two generic invertible operators. By choosing — F2 
S{Q, where S{C) is the squeezing operator 



S{C) = exp 



2„2 



we obtain the tomographic resolution 



(0)= — dxpc{x,(p)TT[OK{X,^c-x)] 

Jo 



(3.47) 



(3.4 



in terms of the probability distribution of the generalized squeezed quadrature 
operators 

X^^ = S^QX^SiC:) = ^ [(^e"^ + ve-''^)a^ + (/ie"*'^ + v*e''P)a\ , (3.49) 

with /i = coshjCI and v = sinh |^| exp[2i arg((^)]. Such an estimation technique 
has been investigated in detail in Ref. |74j. 

A different estimation technique can be obtained by choosing in Eq. H3.46|l 
F\ — I, the identity operator, and F2 — (— )" the parity operator. In this 
case one gets 



O = 



d^a 



Tr 



(3.50) 
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Changing variable to a — 2/3 and using the relation 
it follows 



(3.51) 



(O) = / — Tr 04D\l3){-)'"'D{P) Tr D{f3)pD\f3){~)'"' . (3.52) 



Hence, it is possible to estimate (O) by repeated measurement of the parity 
operator on displaced versions of the state under investigation. An approxi- 
mated implementation of this technique for a single mode radiation field has 
been suggested in Refs. 75» .76, through the measurement of the photon num- 
ber probability on states displaced by a beam splitter. A similar scheme has 
been used for the experimental determination of the motional quantum state of 
a trapped atom |15j . In comparison with the approximated methods, Eq. H3.52|l 
allows to directly obtain the estimator TZ[0]{a) for any operator O for which 
the trace exists. For instance, the reconstruction of the density matrix in the 
Fock representation is obtained by averaging the estimator 



n[\n){n + d\\]ia) = A{n + d\D^ {a){-f ''D{a)\n) 
= 4 (-)"+'^ 



(3.53) 



{n + dy. 



2H^ 



without the need of artificial cut-off in the Fock space ^Hl ■ 



3.3.5 Quantum estimation for spin systems 

The spin tomographic methods of Refs. |2l)[ I69L 1^ allow the reconstruction of 
the quantum state of a spin system. These methods utilize measurements of 
the spin in different directions, i.e. the quorum is the set of operators of the 
form S ■ n, where S is the spin operator and n = (cos(/3sini^, siniySsinT^, cosi?) 
is a varying unit vector. Different quorums can be used, that exploit different 
sets of directions. 

The easiest choice for the set of directions n is to consider all possible di- 
rections. The procedure to derive the tomographic formulas for this quorum is 
analogous to the one employed in Sec. 13.3.31 for homodyne tomography. The 
reconstruction formula for spin tomography for the estimation of an arbitrary 
operator O writes 

(O) = ^ / ^p{m,n)TZ[0]{m,n) , (3.54) 

where p(m, n) is the probability of obtaining the eigenvalue m when measuring 
the spin along direction n, TZ[0]{m, n) is the tomographic estimator for the op- 
erator O, and n is the unit sphere. In this case the operators C(A) of Eq. H3.18|l 



32 



are given by the set of projectors over the eigenstates \m, n) of the operators 
iS* • n. Notice that this is a complete set of operators on the system Hilbert space 
T-L. In order to find the dual basis B, one must consider the unitary operators 
obtained by exponentiating the quorum, i.e. D{ip, ft) — ex.p{iipS ■ n), which sat- 
isfy the bi-orthogonality condition H3.21|l . In fact, D{ip, n) constitutes a unitary 
irreducible representation of the group G = SU(2), and the bi-orthogonality 
condition is just the orthogonality relations between the matrix elements of the 
group representation 1^, i.e. 

[ dg Djr{9)D\^.{g) = ^S^Av , (3.55) 
Jg " 

where Z? is a unitary irreducible representation of dimension d, dg is the group 
Haar invariant measure, and V — J^dg. For G — SU{2), with the 2s -I- 1- 
dimensional unitary irreducible representation D{ip, n) (n £ unit vector on 
the sphere, and tp € [0,47r] the rotation angle around ft) the Haar's invariant 
measure is sin^ ^ sinf} dd dip dip, and ^ — 2s+\ ■ need however to integrate 
only for t/j E [0, 27r] (the change of sign for 27r rotation is irrelevant), whence the 
bi-orthogonality condition writes 

^^^'^ ^ drt r sin^ t 0-|e^^»-^>)(t|e-^'/'"-s>) = 5jk5tr , (3.56) 



47r^ 

and hence the spin tomography identity is given by 



O = 



2s ' 1 '■ '■2'' 



47r2 



- [ dn [ dip sin^ ^ Tr [OD\iP,n)\ D{iP,n) . (3.57) 
Jn Jo 2 



Notice the analogy between Eq. (|3.57|l and Glauber's formula H3.34|l . In fact, 
both homodyne and spin tomography can be derived using the method of Group 
Tomography [201, and the underlying groups are the Weyl-Heisenberg group 
and the SU{2) group, respectively. Formula (|3.54l) is obtained from Eq. H3.57|l 
through the expectation value calculated on the eigenstates oi S ■ ft. Thus, the 
explicit form of the tomographic estimator is obtained as 

n[0]{m, fi) = ^^ii d^ sin2 t Tr \o e-'^M e'^"" . (3.58) 
■K Jo 2 L J 

As already noticed, there are other possible quorums for spin tomography. 
For example, for spin s = ^ systems, a self-dual basis for the operator space 
is given by the identity and the Pauli matrices. In fact, from the proper- 
ties Tr[(TQ,] = and CTqCT/j = ij^j ^ajS^yO'aO'p {a, (3,^ = x,y,z), both the bi- 
orthogonality relation (|3.21() and the trace condition H3.23|l follow. In this case 
the reconstruction formula writes 

(O) - ^Tr [O] + i ^ Tr [Oa,] . (3.59) 
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In the case of generic s spin system, Weigert has also shown ^\ that by 
choosing (2s + 1)^ arbitrary directions for n, it is possible to obtain (in almost 
all cases) a quorum of projectors \s, nj){s, nj\ {j — 1, - ■ ■ , (2s + l)^), where \s, fij) 
is the eigenstatc pertaining to the maximum eigenvalue s of 5 • fij . 

3.3.6 Quantum estimation for a free particle 

The state of a moving packet can be inferred from position measurement at 
different times |7S]. Assuming a particle with unit mass and using normalized 
unit h/2 = 1, the free Hamiltonian is given by the square of momentum operator 
Hp = ■ In terms of the eigenvectors |a;) of the position operator and of the 
selfadjoint operator 

R{x,t) = e-'P^^ \x){x\ e'f'^ , (3.60) 

the probability density of the position of the free particle at time r is obtained 
as p(x,t) — Tt:[p R[x,t)]. The operators R{x,t) provide a self-dual basis, and 
an arbitrary particle state can be written as 

P = I I dx dr p{x,t) R{x,t) . (3.61) 
Js. Jr 

3.4 Noise deconvolution and adaptive tomogra- 
phy 

In this section we will analyze: 1) the noise deconvolution scheme of Refs. |2()l 
I79j . that allows to eliminate the experimental noise that arises from imperfect 
detection and lossy devices; 2) the adaptive tomography technique of Ref. [22] 
that allows to tune the unbiased tomographic estimators to a specific sample of 
experimental data, in order to reduce the statistical noise. 

3.4.1 Noise deconvolution 

In short, it is possible to eliminate detection noise when it is possible to invert 
the noise map. A noise process is described by a trace preserving completely 
positive map T. The noise can be deconvolved at the data analysis if 

• the inverse of T exists, namely T^^ : /:(H) C{n), with T^^ [T [O]] = O, 

for yoecin). 

• the estimator £\[0]{Q\) is in the domain of F^^ 

• the map T^^ [£x[0]{Qx)] is a function of Q\. 

If the above conditions are met, we can recover the "ideal" expectation value 
(O) that we would get without noise. This is achieved by replacing £\[0]{Q\) 
with r~^[£A[0]((3A)], and evaluating the ensemble average with the state r'^(p). 
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namely the state affected by the noise {V^ represents the dual map, that provides 
the evolution in the Schroedinger picture). Hence, one has 

(O)- / dXTT[r-^[£),[0]{QxW{p)]= / dA(r-i[^A[0](gA)])r . (3.62) 



Consider for example the noise arising from non-unity quantum efficiency r] 
of homodyne detectors. Recall that the ideal probability density is replaced by 
a Gaussian convolution with rms A^^ = (1 — 77)/(4?7). Then, the map F,, acts on 
the quorum as follows 



f + OO 

r^[e^*=^^] = / dxe'''^r^[\x){x\] (3.63) 



— OO 
+ 00 /' + 00 



/ ikx - \\^f\ /_/n -hA^k^ ikX, 



dx / dx'e'^^'e — ^[\x'){x'\]^e-^ 

J —oo J~OQ 

Of course one has 

^giA^fe^g^fcx^ ^ (3.64) 

In terms of the Fourier transform of the estimator 

r+°° dx 

n[0]{y,ip)^J — e"^7^[0](x,^p) , (3.65) 

one has 

TZMy, = e^^'^'7^[0](y, ^) . (3.66) 

We applied the above result in Sec. 13. 3. 31 where the effect of non-unity quantum 
efficiency for reconstructing the density matrix elements was discussed. The use 
of the estimator in Eq. H3.42(l and the origin of the bound rj > 1/2 is now more 
clear. 

Another simple example of noise deconvolution is given here for a spin 1/2 
system. Consider the map that describes the "depolarizing channel" 

Tp[0]^il-p)0 + ^TT[0]1 , 0<p<l. (3.67) 
This map can be inverted for p 7^ 1 as follows 

r,-MO] = Y^(0-|Tr[0]/) . (3.68) 
Then Eq. (|3.59|) can be replaced with 

(O) = ^Tr [O] + E "^Pf ("^' , (3.69) 

m=±^ a=x,y,z 

where now Pp{m, n^) represents the probability of outcome m when measuring 
aa on the noisy state Fp[p]. 
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3.4.2 Adaptive tomography 

The idea of adaptive tomography is that the tomographic null estimators of Eq. 
H3.30|l can be used to reduce statistical errors. In fact, the addition of a null 
estimator in the ideal case of infinite statistics does not change the average since 
its mean value is zero, but can change the variance. Thus, one can look for a 
procedure to reduce the variance by adding suitable null functions. Consider 
the class of equivalent estimators for O 

M 

S'x[0]{Qx) = £x[0]{Qx) + J2 ''^MQx) ■ (3.70) 

i=l 

Each estimator in the class £' is identified by the coefficient vector i7. The 
variance of the tomographic averages can be evaluated as 

M M 

A2£'[0] = A2£[0] + 2 ^ l^^M,£[0] + J2 > (3-71) 

i=l i,j=l 

where F= {Jj^dX F{Qx)), and 

Minimizing A2£'[0] with respect to the coefficients I'i, one obtains the equation 

M 

^j^^j = -WW^ , (3.73) 

which can be solved starting from the estimated mean values, with the vector v 
as unknown. Notice that the obtained vector v will depend on the experimental 
data, and has to be calculated with the above procedure for any new set of data. 

In this way we obtain an adaptive tomographic algorithm, which consists in 
the following steps: 

• Find the null estimators Mi{Qx) (* ~ li ■ ■ ■ 7 M) for the quorum which is 
being used in the experiment. 

• Execute the experiment and collect the input data. 

• Calculate, using the obtained data, the mean values AfiAfj and £[0]Afi, 
and solve the linear system (|3.73|) . to obtain i7. 

• Use the vector v obtained in the previous step to build the 'optimized 
estimator' £'[0]{Qx) = £[0]{Qx)+J2z ^i^i{Q\)- Using the data collected 
in the first step, the mean value (O) is now evaluated as 

(O) = [ d\ {E'x[0]{Qx)) , (3.74) 

where the optimized estimator has been used. 
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• For each new set of data the whole procedure must be repeated, as v is 
dependent on the data. 

Notice that also the experimental mean values are slightly modified in the adap- 
tive tomographic process, since null estimators do not change mean values only 
in the limiting case of infinite statistics. Examples of simulations of the adap- 
tive technique that efficiently reduce statistical noise of homodync tomographic 
reconstructions can be found in Ref. |22j . In homodync tomography null esti- 
mators are obtained as linear combinations of the following functions 

AfkAX^) = K e±'('=+2+'")'^ , fc, n > . (3.75) 

One can easily check that such functions have zero average over (p, indepen- 
dently on p. Hence, for every operator O one actually has an equivalence class 
of infinitely many unbiased estimators, which differ by a linear combination 
of functions N'kjuiX^). It is then possible to minimize the rms error in the 
equivalence class by the least-squares method, obtaining in this way an optimal 
estimator that is adapted to the particular set of experimental data. 
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Chapter 4 

Universal homodyning 



As shown in Ref. J^, homodyne tomography can be used as a kind of univer- 
sal detector, for measuring generic field operators, however at expense of some 
additional noise. In this chapter, the general class of field operators that can 
be measured in this way is reviewed, which includes also operators that are 
inaccessible to heterodyne detection. In Ref. the most relevant observ- 

ables were analyzed — such as the intensity, the real, the complex field, and the 
phase — showing how their tomographic measurements are affected by noise that 
is always larger than the intrinsic noise of the direct detection of the considered 
observables. On the other hand, by comparing the noise from homodyne tomog- 
raphy with that from heterodyning (for those operators that can be measured 
in both ways), in Ref. 29 it was shown that for some operators homodyning is 
better than heterodyning when the mean photon number is sufficiently small, 
i.e. in the quantum regime, and in this chapter such comparison will be also 
reviewed. 



4.1 Homodyning observables 

Homodyne tomography provides the maximum achievable information on the 
quantum state of a single-mode radiation field through the use of the estimators 
in Sec. 13.3.31 In principle, the knowledge of the density matrix should allow one 
to calculate the expectation value also for unbounded operators. However, this 
is generally true only when one has an analytic knowledge of the density matrix, 
but it is not true when the matrix has been obtained experimentally. In fact, 
the Hilbert space is actually infinite dimensional, whereas experimentally one 
can achieve only a finite matrix, each element being affected by an experimental 
error. Notice that, even though the method allows one to extract any matrix 
element in the Hilbert space from the same bunch of experimental data, it is the 
way in which errors converge in the Hilbert space that determines the actual 
possibility of estimating the trace (O) ~ Tr[Op] for an arbitrary operator O. 
This issue has been debated in the set of papers of Ref. |73]- Consider for 
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example the number representation, and suppose that we want to estimate the 
average photon number {a''' a). In Ref. |Hni it has been shown that for non- 
unit quantum efficiency the statistical error for the diagonal matrix element 
{n\p\n) diverges faster than exponentially versus n, whereas for rj — I the error 
saturates for large n to the universal value £„ = •y/2/iV that depends only on 
the number N of experimental data, but is independent on both n and on the 
quantum state. Even for the unrealistic case 77 = 1, one can see immediately that 
the estimated expectation value {a' a) = "^^=0 ^^Pnn based on the measured 
matrix elements p„„, will exhibit an unbounded error versus the truncated- 
space dimension H, because the non- vanishing error of p„„ versus n multiplies 
the increasing eigenvalue n. 

Here, we report the estimators valid for any operator that admits a normal 
ordered expansion, giving the general class of operators that can be measured in 
this way, also as a function of the quantum efficiency 77. Hence, from the same 
tomographic experiment, one can obtain not only the density matrix, but also 
the expectation value of various field operators, also unbounded, and including 
some operators that are inaccessible to heterodyne detection. However, the 
price to pay for such detection flexibility is that all measured quantities will be 
affected by noise. If one compares this noise with that from heterodyning (for 
those operators that can be measured in both ways), it turns out that for some 
operators homodyning is anyway less noisy than heterodyning, at least for small 
mean photon numbers. The procedure for estimating the expectation (O) will 
be referred to as homodyning the observable O. 

By homodyning the observable O we mean averaging an appropriate esti- 
mator TZ[0]{x, If), independent on the state p, over the experimental homodyne 
data, achieving in this way the expectation value (O) for every state p, as in 
Eq. (|3.3t)|) . For unbounded operators one can obtain the explicit form of the 
estimator TZ[0]{x, if) in a different way. Starting from the identity involving 
trilinear products of Hermite polynomials j81| 

dxe Hk(x)Hm{x)Hn{x) = —- - , (4.1) 

-00 {s - k)\{s - m)\{s - ny. 

for fc + m + rt = 2s even, Richter proved the following nontrivial formula for the 
expectation value of the normally ordered field operators 



„tn„™^^/ '^Jil rf^p(^,^)e-("'-»)^-"+"V-^-;^^^ , (4.2) 



which corresponds to the estimator 

7^[at«a™](.T,(^) = . (4.3) 

This result can be easily extended to the case of non-unit quantum efficiency 
r] <l. Using Eq. I|3.(j6|l one obtains 



,^7 ^(2,y)n+™(«+™) 
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From Eq. H4.4|l by linearity one can obtain the estimator TZri[f]{x, (p) for any 
operator function / that has normal ordered expansion 



nm— 



From Eq. H4.4|) one obtains 



(N) in m 
nm " " 



(4.5) 



+m,s 



E 

s=0 



i;=0 



where 



nm— 

Continuing from Eq. (|4.()|l one has 



(4.6) 



(4.7) 



7^»)[/](a;,</5) = exp 



2?] du^ 



212: d 
y/fj dv 



v=0 



and finally 



^»?[/](2;,</') 



dw 



T[,f]{w + 2ix/^, ^) 



(4.8) 



(4.9) 



Hence one concludes that the operator / can be measured by homodyne tomog- 
raphy if the function jF[/](v, ip) in Eq. (|4.7|) grows slower than exp(— 7711^/2) for 
D — » 00, and the integral in Eq. I|4.9|) grows at most exponentially for x ^ 00 
(assuming p(x, ^p) goes to zero faster than exponentially at x — + 00). 

The robustness to additive phase-insensitive noise of this method of homo- 
dyning observables has also been analyzed in Ref. |16 | , where it was shown that 
just half photon of thermal noise would spoil completely the measurement of 
the density matrix elements in the Fock representation. 

In Table UTI we report the estimator 7?,^[0](a;, (/?) for some operators O. 
The operator Ws gives the generalized Wigner fimction M4(a, a*) for order- 
ing parameter s through the relation in Eq. I|2.11(l . From the expression of 
7?,^[Ws](a;, if) it follows that by homodyning with quantum efficiency 77 one can 
measure the generalized Wigner function only for s < 1 — 77^^: in particular the 
usual Wigner function for s = cannot be measured for any quantum efficiency. 
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(ata)2 



7r(l-s) \s-l 



dt 



2e- 



^(1 - - ^ 



cos 2, 



2t 



7^^[|n)(n + c^|](a;,y>) in Eg. (|3.44l 



jn) (n + d\ 



Table 4.1: Estimator 7^,,[0](x, (/?) for some operators O (From Ref. ll6| 1. 



4.2 Noise in tomographic measurements 

In this section we will review the analysis of Ref . , where the tomographic 
measurement of following four relevant field quantities have been studied: the 
field intensity, the real field or quadrature, the complex field, and the phase. 
For all these quantities the conditions given after Eq. (|4.9|l are fulfilled. 

The tomographic measurement of the observable O is provided in terms of 
the average of the estimator = TZri[0]{x, if) over the homodyne data. The 

precision of the measurement is given by the confidence interval 
is a real quantity, one has 



Aii;2 = u>2 -iZy^2 ^ (4.10) 

where 

^ = ^2p](^^ /■ dxp^{x,ip)nl[0]{x,ip) . (4.11) 

Jo ^ J-co 

When is complex, one has to consider the eigenvalues of the covariance 
matrix, namely 



1 



(4.12) 

When the observable O can also be directly measured by a specific setup we 
can compare the tomographic precision Aw^ with (AO^) — (O^) — (O)^. Notice 
that, when we deal with rj < 1 the noise (AO^) is larger that the quantum 
fluctuations due to smearing effect of non-unit quantum efficiency. As we will 
see, the tomographic measurement is always more noisy than the corresponding 
direct measurement for any observable at any quantum efficiency rj. This is 
not surprising, in view of the larger amount of information retrieved in the 



41 



tomographic measurement as compared to the direct measurement of a single 
quantity. 

According to Eq. (|4.11|) . the evaluation of the added noise requires the 
average of the squared estimator. For the estimators in Eq. 14. 4|) it is very 
useful the following identity for the Hermite polynomials [53] 

that allows to write 

7^2[a^"a™](a;,(^) = 

^ y ^M a^^'Ai^M , (4.14) 



fc=0 



namely the squared estimator 7?.^[a^"a'"](x, (y9) can be written just in terms of 
"diagonal" estimators 'Kj^\a}^ a^\{x , ip). 

4.2.1 Field-Intensity 

Photodetection is the direct measurement of the field-intensity. For non-unit 
quantum efhciency ij, the probability of detecting m photons is given by the 
Bernoulli convolution in Eq. (|2.22() . Let us consider the rescaled photocurrent 

= ifl^a , (4.15) 

which traces the photon number, namely 

(/,,) — — m Priini) — {a' a) = n . (4-16) 

The variance of is given by 

1 °° /I \ 

(A/2) - — 51 ™^P(w) - = (An^) + fi - - 1 , (4.17) 



' m=0 ^ ' 



where (An^) denotes the intrinsic photon number variance, and n{rj~^ — 1) 
represents the noise introduced by inefficient detection. The tomographic esti- 
mator that traces the photon number is given by the phase-independent function 
Wri = 2x^ — {2r])~^. Using Eq. (|4.14ll we can evaluate its variance as follows 

1,9, _ /2 3\ 1 
" 2 j + V 



Aw?, = {An') + -{v?) +n[l-"-)+—. (4.18) 



The noise iV[n] added by tomography in the measurement of the field intensity 
n is then given by 



N[n] = Awl - = I 







1 " 


in') 













(4.19) 
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Notice that N[n] is always positive, and largely depends on state under exami- 
nation. For coherent states we have the noise-ratio 



Srirj — 








1/2 










rjn J 





(4.20) 



4.2.2 Real Field 

For single mode radiation the electric field is proportional to a quadrature X = 
(a + a^)/2, which is just traced by homodyne detection at fixed zero-phase with 
respect to the local oscillator. The tomographic estimator is given by Wj) = 
TZri[X]{x, (p) = 2a; cos (/9, independently on t], whereas the squared estimator 
7?.^[X] can be written as 

= 1 [7^^ [a^] (x, ^) + [a^^] (x, ip)] + TZ^ [a^a] (x, v) + ^ + ^^^^ • (4.21) 



Then 



has 



Aw2 



1 

2^ 



(AX' 



(4.22) 



where (AX'^) represents the intrinsic quadrature fluctuations. The tomographic 
noise in Eq. H4.22(l can be compared with the rms variance of direct homodyne 
detection (see Sec. I2.3|l 
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(AX'), = (AX^) 
Then the added noise reads 

r , n 1 
^[^]-2 + 4^- 

For coherent states (AX^) = 1/4, and one has the noise-ratio 



(4.23) 



(4.24) 



Sx,j — 



rjn 



(4.25) 



4.2.3 Field amplitude 

The detection of the complex field amplitude of a single-mode light-beam is rep- 
resented by the generalized measurement of the annihilation operator a. The to- 
mographic estimator for a is given by the complex function = 7?.^[a](x, (p) = 
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2xexp(iip), and the precision of the measurement is evaluated as in Eq. (I4.12|) . 
From Eq. 14.14|l one obtains 



= 7^^[a](x,^p) 

and 



— h 2TZ,j[a^ a]{x , f ) 



— +n,j[a^]ix,ip) , (4.26) 



kijl^ = |7^»?H(a;, <^)|^ = - [1 + 2777^„[a^a](x, , 



(4.27) 



and hence 



- + 2n~\(a)\^±\{a^) 



(4.28) 



The optimal measurement of the complex field a is obtained through heterodyne 
detection. As notice in Sec. 12.41 the probability distribution is given by the 
generalized Wigner function Ws{a,a*), with s = 1 — ^. Using Eq. H2.56|l the 
precision of the measurement is easily evaluated as follows 



— 



n + --\{a)\'±\{a')~{ay 



V 



(4.29) 



The noise added by quantum tomography then reads 

N[a] = in , 



(4.30) 



which is independent on quantum efficiency. For a coherent state we have 



Aw;2 = i 



1 

n -\ — 
V. 



1 



and the noise ratio then writes 
6a„ = \ 



<Aa% 



\/l + r]n 



(4.31) 



(4.32) 



4.2.4 Phase 

The canonical description of the quantum optical phase is given by the proba- 
bility operator measure |84l 



dip 
2^ 



exp[z(TO — n)ip] \n) {m\ 



(4.33) 



n,m— 



However, no feasible setup is known that achieves the optimal measurement 
H4.33|l . For this reason, here we consider the heterodyne measurement of the 
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phase, and compare it with the phase of the tomographic estimator for the 
corresponding field operator a, i.e. Wjj = arg(2a;e*'''). Notice that the phase Wr, 
does not coincide with the local oscillator phase ip, because x has varying sign. 
The probability distribution of can be obtained by the following identity 



^ dip 



which implies 



dx p-qix, ip) = 1 = 



^ dwn 



dx Pri{x, Wri) 



dx p-q{x, Wri) 



(4.34) 



(4.35) 



The precision in the tomographic phase measurement is given by the rms vari- 
ance Aw^ of the probability (|4.35|) . In the case of a coherent state with positive 
amplitude = ||/3|), Eq. gives 



Pvi^v) 



1 

2^ 



1 + Erf 



■ 721/31 



cos Wn 



Vv 



(4.36) 



which approaches a "boxed" distribution in [—tt/2, 7r/2] for large intensity \(3\ ^ 
1. We compare the tomographic phase measurement with heterodyne detection, 
namely the phase of the direct-detected complex field a. The outcome probabil- 
ity distribution is the marginal distribution of the generalized Wigner function 
W,ia,a*){s = l 



jj) integrated over the radius 



p„(^)= / pdpWs{pe''^,pe-'n 
Jo 



(4.37) 



whereas the precision in the phase measurement is given by its rms variance 
We are not able to give a closed formula for the added noise N[(p] = 

= ||/3|) (zero mean 

7rVl2 and 



— Aip'^. However, for high excited coherent states 



phase) one has Awfj = 
is thus given by 

Siprj 



(27771) . The asymptotic noise-ratio 



A^2 



(4.38) 



A comparison for low excited coherent states can be performed numerically. 
The noise ratio Sip^j (expressed in dB) is shown in Fig. 14. II for some values of 
the quantum efficiency rj. 

It is apparent that the tomographic determination of the phase is more noisy 
than heterodyning also in this low-intensity regime. 

In Table a synthesis of the results of this section is reported. We have 
considered the ratio between the tomographic and the direct-measurement noise. 
This is an increasing function of the mean photon number n, however scaled by 
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5 



Figure 4.1: Ratio between tomographic and heterodyne noise in the measure- 
ment of the phase for fow excited coherent states, The noise ratio is reported 
versus the mean photon number n for some values of the quantum efficiency. 
From bottom to top we have 77 = 0.2, 0.4, 0.6, 0.8, 1.0 (From Ref. fl^). 
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Table 4.2: Added noise N[0] in tomographic measurement of O and noise ratio 
(50,,for coherent states. For the phase the results are valid in the asymptotic 
regime 1 (From Ref. [211 )• 

the quantum efhciency rj. Therefore homodyne tomography turns out to be a 
very robust detection scheme for low quantum efficiency. 

In Fig. 14.21 the coherent-state noise ratio (in dB) for all the considered 
quantities are plotted for unit quantum efficiency versus h. 

In conclusion, homodyne tomography adds larger noise for highly excited 
states, however, it is not too noisy in the quantum regime of low h. It is then 
very useful in this regime, where currently available photodetectors suffer most 
limitations. Indeed, it has been adopted in recent experiments of photodetection 

[iniin]. 
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5 10 

5 

Figure 4.2: The coherent-state noise ratio (in dB) for all the quantities consid- 
ered in this section (From Ref. j29j). 

4.3 Comparison between homodyne tomography 
and heterodyning 

We have seen that homodyne tomography allows to measure any field ob- 
servable / = f{a,a^) having normal ordered expansion / = /^^^(a, a^) = 

Zl^=o Z"™'*'^^""'" ^^"^ bounded integral in Eq. (|4.9|l . On the other hand, 
as shown in Sec. 12.41 heterodyne detection allows to measure field observables 
that admit anti-normal ordered expansion f = J^^\a, a)) = X]^=o fvmCL^a^^, 
in which case the expectation value is obtained through the heterodyne average 

(/)- / —f^^\a,a*){a\p\a) . (4.39) 

As shown in Sec. 12.41 for 77 = 1 the heterodyne probability is just the Q- 
function Q(a, a*) — ^(alplQ;), whereas for 77 < 1 it is Gaussian convoluted with 
rms (1 — Ty)/??, thus giving the Wigner function Wa(a^ a*), with s = 1 — |. 

Indeed, the problem of measurability of the observable / through heterodyne 
detection is not trivial, since one needs the admissibility of anti-normal ordered 
expansion and the convergence of the integral in Eq. I|4.39|l . We refer the reader 
to Refs. ^SF, 'TF for more details and to Refs. [511 ED) for analysis of quantum 
state estimates based on heterodyne detection. 

The additional noise in homodyning the complex field a has been evaluated 
in Eq. (|4.3U|) . where we found that homodyning is always more noisy than 
heterodyning. On the other hand, for other field observables it may happen that 
homodyne tomography is less noisy than heterodyne detection. For example, 
the added noise in homodyning the intensity a with respect to direct detection 
has been evaluated in Eq. (|4.19() . Analogously, one can easily evaluate the added 
noise ]>lhf,t\n\ when heterodyning the photon number n — a. According to Eq. 
I|2.5()|l . the random variable corresponding to the photon number for heterodyne 
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detection with quantum efficiency 77 is i^(a) = \a\^ — ^. From the relation 

\^ = (a2at2) + (aat) + 2 f (4.40) 

V \ V J 



one obtains 

= (An2) + n f - - 1) + J- . (4.41) 

Upon comparing with Eq. H4.17|) , one concludes that the added noise in hetero- 
dyning the photon number is given by 



Nhet[n] = Aiy^iz) - {Al') ^—{r]n-l). (4.42) 
With respect to the added noise in homodyning of Eq. 14.19|l one has 

Nhet[n] = N[n] - ^ ( (n^) - n - ^ ) . (4.43) 



2 V V 

Since (n^) > we can conclude that homodyning the photon number is 
less noisy than heterodyning it for sufficiently low mean photon number (n) < 
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Chapter 5 



Multimode homodyne 
tomography 

The generalization of homodyne tomography from a single-mode to a multimode 
field is quite obvious, the estimator of simple operator tensors O ^ Oi ® O2 ® 
. . .® On being just the product of the estimators of each single-mode operator 
Oi , Oi , . . . , On ■ By linearity, one then obtains also the estimator for arbitrary 
multimode operators. Such a simple generalization, however, requires a separate 
homodyne detector for each mode, which is unfeasible when the modes of the 
field are not spatio-temporally separated. This is the case, for example of pulsed 
fields, for which a general multimode tomographic method is especially needed, 
also due to the problem of mode matching between the local oscillator and the 
detected fields (determined by their relative spatio-temporal overlap) |H3|) which 
produces a dramatic reduction of the overall quantum efficiency. 

In this Chapter we review the general method of Ref. jl7| for homodyning 
observables of a multimode electromagnetic field using a single local oscillator 
(LO), providing the rule to evaluate the estimator of an arbitrary multimode op- 
erator. The expectation value of the operator can then be obtained by averaging 
the estimator over the homodyne outcomes that are collected using a single LO 
whose mode randomly scans all possible linear combinations of incident modes. 
We will then specifically consider some observables for a two-mode field in a 
state corresponding to a twin-beam produced by parametric downconversion, 
and prove the reliability of the method on the basis of computer simulations. 

Finally, we report some experimental results Efil obtained in the Prem Ku- 
mar's lab at Northwestern University. Such experiment actually represents the 
first measurement of the joint photon-number probability distribution of the 
twin-beam state. 
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5.1 The general method 

The Hilbert-Schmidt operator expansion in Eq. H3.34|l can be generalized to 
any number of modes as follows 



O = 



(P Zq f Cp Zl 



(PzM 



Tr <^ O exp 



■ M 
.1=0 



Zl ai 



X exp 



■ M 
./=0 



Zl ai 



(5.1) 



where a; and aj, with I = 0, . . . ,M and [a;, aj,] = Sw , are the annihilation and 
creation operators of AI +1 independent modes, and O now denotes an operator 
over all modes. Using the following hyper-spherical parameterization for € C 



(5.2) 



zo='-kuQ{9)e'^'> = 


^fce^'^o cos 6*1 , 


z,^'-ku,{9)e'^' = 


^fce"^i sin6li cos 6*2 , 




^ k e^^^ sin 9i sin 92 cos 


ZM-i='i^kuM-i{e)e'^^'-^ = 


^fce*'^'^-i sin6'i sin 6*2. 


ZM^\kuM{e)e'^'' = 


Uie''^^' sin6li sin6l2 . . 


where k e [0,oo); V'z e [0, 27r] for I = 0,1,..., M; and 
1,2,..., M, Eq. IjS.lfl can be rewritten as follows: 





Here we have used the notation 



+00 /,x2M+l . 

\ 2 / M! ^ 



dipi 



df^m = n 

y" d^[9] = 2^' M\ |] fj^ d9i sin2(^^-')+i 0, cos 0, 
X{9,{,) = \[a\1^) + A{9,^) 

M 

A(0;V^) = ^e-^'^'ziz(%z. 



(5.3) 

(5.4) 

(5.5) 
(5.6) 
(5.7) 



1=0 



From the parameterization in Eq. (|5.3(l . one has X]i!lo"?(^) ~ 1; hence 
[A(0, ^/;), >1^"(6', -0)] = 1, namely A{9, tp) and A^{9, -0) themselves are annihilation 
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and creation operators of a bosonic mode. By scanning all values of 9i G [0, it/2] 
and "0; G [0, 27r], all possible linear combinations of modes ai are obtained. 

For the quadrature operator X{9,ip) in Eq. ()5.6|l . one has the following 
identity for the moments generating function 

{e^'^''^''^'') = exp (^fc') f°°dxe^''-p,{x;ej) , (5.8) 

where Pri{x; 9, ip) denotes the homodyne probability distribution of the quadra- 
ture X{9,'tjj) with quantum efhciency r/. Generally, r/ can depend on the mode 
itself, i.e., it is a function 77 ~ 'ri{9,ip) of the selected mode. In the following, 
for simplicity, we assume t] to be mode independent, however. By taking the 
ensemble average on each side of Eq. (15. 3|) and using Eq. (|5.8|l one has 

r r r+co 

{O) = d^i[i^] dfi[9] / dxp,j{x;9,i^)nrj[0]{x;9,^) , (5.9) 



where the estimator TZ.,j[0]{x; 9, ijj) has the following expression 

n„[0]{x;9j) ^ —— / dte-(i-t)*+2»v^-t*^Tr[Oe-2'^^(^''^)],(5.10) 
^I- Jo 

with K = 2-q/ {2-q—l). Eqs. H5.9|l and H5.10|l allow to obtain the expectation value 
(O) for any unknown state of the radiation field by averaging over the homodyne 
outcomes of the quadrature X{9, ip) for 9 and -0 randomly distributed according 
to diJ,[ip] and dfi[9]. Such outcomes can be obtained by using a single LO that is 
prepared in the multimodc coherent state <Xif£olT') '^itli 7/ — e^'^'ui{9)K/2 and 
K ^ 1. In fact, in this case the rescaled zero-frequency photocurrent at the 
output of a balanced homodyne detector is given by 

1 

(=0 

which corresponds to the operator X{9, In the limit of a strong LO {K 
00), all moments of the current / correspond to the moments of X{9,ip), and 
the exact measurement of X{9,il;) is then realized. Notice that for modes ai 
with different frequencies, in the d.c. photocurrent in Eq. H5.11(l each LO with 
amplitude 7; selects the mode ai at the same frequency (and polarization). For 
less-than-unity quantum efficiency, Eq. (|5.8|) holds. 

Equation H5.10|l can be applied to some observables of interest. In particular, 
one can estimate the matrix element ({ni}|_R|{mi}) of the multimodc density 
operator R. This will be obtained by averaging the estimator 

'Jl,[\{mi}){{ni}\]{x;9,i;) = e-S,.o(n,-™,)V, __ 
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r-\-oo ^'^ 

X dte~*+2*^^i^'^+S"o(M<-^i)/2-Q2,;;i-''<[Ku2(0)t] , (5.12) 

■'O 1=0 

where = max(mi, n;), i^i = min(mi, n;), and L'^{z) denotes the generalized 
Laguerre polynomial. For diagonal matrix elements, Eq. (|5.12|l simplifies to 



' z=o 

with Ln{z) denoting the customary Laguerre polynomial in z. Using the follow- 
ing identity [HI] 

^ao+ai+...+a„+Af + Xl + . . . + Xm) 

E ir(^o)ir/(xi)...Lf-(xM), (5.14) 

io+il + ---+iM=n 

from Eq. H5.13|l one can easily derive the estimator of the probability distribution 
of the total number of photons N — X]i=o '^I'^i 

n.,[\n){n\]{x;0,^) = dt e-'+'^^' ^ t'' L^' [>,t] , (5.15) 

where |n) denotes the eigenvector of N with eigenvalue n. Notice that the 
estimator in Eq. (|5.13|l does not depend on the phases ■(/;;; only the knowledge 
of the angles 9i is needed. For the estimator in Eq. (|5.15(l . even the angles 9i 
can be unknown. 

Now we specialize to the case of only two modes a and b (i.e., M=l and 9 
is a scalar 9). The joint photon- number probability distribution is obtained by 
averaging 



TZrj[\n,m){n,m\]{x;9,'>po,i^i) = 

r+oo 

dte-^^+'^'^^'' t Lnint cos^ 9)L^{Ktsw? 9) . (5.16) 



The estimator (|5.15|l of the probability distribution of the total number of pho- 
tons can be written as 

nr,[\n){n\]{x;9,^o,tpi) ^ dt e'*+^"^' t L^lnt] . (5.17) 

Jo 

For the total number of photons one can also derive the estimator of the moment 
generating function, using the generating function for the Laguerre polynomials 
|81|. One obtains 

a+a+b+hi/ . X _ 1 ^ / „ 1. ^~ z 



7^,[z"''^+'''''](x;0,Vo>l) = — -T^* 2,-; -j^ x' ] . (5.18) 
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For the first two moments one obtains the simple expressions 

7^,Ja^a + 6^6](x;6',V'o,V'l) ='ix^ + --2, (5.19) 

nj{a)a + b'^bf]{x- d, V-o, V'l) = 8x^ + (— - 2o) + ~ — + 4 . 

V 7 / 7 7 

It is worth noting that analogous estimators of the photon-number difference 
between the two modes are singular and one needs a cutoff procedure, similar 
to the one used in Ref. [SZ] for recovering the correlation between the modes 
by means of the customary two- mode tomography. In fact, in order to extract 
information pertaining to a single mode only one needs a delta- function at 9 = 
for mode a, or 9 = tt/2 for mode b, and, in this case, one could better use 
the standard one-mode tomography by setting the LO to the proper mode of 
interest. 

Finally, we note that for two-mode tomography the estimators can be aver- 
aged by the integral 

{O) = -7^-1 -7^-1 ^ / dxp^{x;9,il)o,i}i) 

Jq Jo J-1 2 

X 7^,,[0] (a:; 0,^0^1) (5.20) 

over the random parameters cos(20), ^/jq, and V'l- For example, in the case of 
two radiation modes having the same frequency but orthogonal polarizations, 
9 represents a random rotation of the polarizations, whereas "00 and "01 denote 
the relative phases between the LO and the two modes, respectively. 



5.1.1 Numerical results for two- mode fields 

In this section we report some Monte-Carlo simulations from Ref. Tf to judge 
the experimental working conditions for performing the single-LO tomography 
on two-mode fields. We focus our attention on the twin-beam state, usually 
generated by spontaneous parametric downconversion, namely 

oo 

1*) = 5(x)|0)a|0), = v/T3^^^" \n)a\n), , (5.21) 

n=0 

where S{x) — exp(xa^6''' — x*o-b) and ^ = e*'*''^-^ tanh|x|. The parameter ^ 
is related to the average number of photons per beam n = W\^/{^ — ICP)- 
For the simulations we need to derive the homodyne probability distribution 
p{x; 9, 'ipQ^ipi) which is given by 

pix;9,^Po,i'i) = Tr[UUx)aa{x\®lbU\^)m] (5.22) 

= ,(0|fc(0| S^X) [\x)aa{x\ <E> U]US{X) \0)a\0)b , 

where \x)a is the eigenvector of the quadrature x = ^{a^ + a) with eigenvalue x 
and U is the unitary operator achieving the mode transformation 

^ [b)^=[-e^^^sm9 e^^»cos9 ) [b) ' ^^'^^^ 
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In the case of two radiation modes having the same frequency but orthogo- 
nal polarizations — the case of Type II phase-matched parametric amplifier — 
Eq. H5.22|l gives the theoretical probability of outcome x for the homodyne 
measurement at a polarization angle 9 with respect to the polarization of the 
a mode, and with ipQ and ■01 denoting the relative phases between the LO 
and the two modes, respectively. By using the Dirac-(5 representation of the 
X-quadrature projector 



-exp[zA(X-x)], 



(5.24) 



Eq. (|5.22(l can be rewritten as follows ^Jj 



p(x;6',V'o, V'l) 



27r 



iXx 



27r 



I exp < I 



.A 



^e--'V* cos 61 -t- e~*'''Vsine')6 + H.c.]} |0)a|0)b 
where we have used Eq. H5.23|l and the transformation 



fj, V 



with \x = cosh|x| and v = e*"8;x sinh|x|. Upon defining 



KC ^ e-''^''^ cost 



KD^e'^^^v* cosO + e ->smf , 
where € R and C, L» G C, with [Cp + |D|2 ^ l one has 

if^ = + + 2fi\iy\ sin26'cos(V'o + V'l - argj/) 
Now, since the unitary transformation 



C D 

D* C* 



(5.25) 
(5.26) 

(5.27) 
(5.28) 

(5.29) 



has no effect on the vacuum state, Eq. (|5.25|) leads to the following Gaussian 
distribution 



p(a;;6','0o,'0i) 

2^ 
d\ 
2^ 



—iXi 



— oo 

+ 00 



|exp<jzif- [iCa + Db) + Il.c.]}\0)a\0)b 



V27rA2(0,V'o,V'i 



I exp 
: exp 



iK- (a + a^) 



|0)a- 



K 



<^/K)a\ 



2A2(0,V'o,^i) 



(5.30) 
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where the variance A^(0, -00, V"!) is given by 



_ l + |^|2 + 2|^|sin26lcos(V'o + V'i-argO 



4(1 -ICP) 



(5.31) 



Taking into account the Gaussian convolution that results from less-than-unity 
quantum efficiency, the variance just increases as 



A2(0,Vo>i) ^ A'(^,0o,0i) = A'(0,0o,V^i) + 



1 — 77 
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(5.32) 



Notice that the probability distribution in Eq 
vacuum for = j and -00 + "01 ^ arg = or tt 



corresponds to a squeezed 



p (n, m) 




p (n, m) 




-0 . 05 



Figure 5.1: Two-mode photon-number probability p{n,m) of the twin-beam 
state in Eq. H5.21|l for average number of photons per beam n = 5 obtained by a 
Monte-Carlo simulation with the estimator in Eq. (|5.16ll and random parameters 
cos 20, ipo, and ■01- On the left: quantum efficiency 77 = 1 and 10^ data samples 
were used in the reconstruction. On the right: 77 — 0.9, and 5 x 10^ data samples 
(From Ref. [TT]). 

We study the tomographic measurement of the joint photon-number proba- 
bility distribution and the probability distribution for the total number of pho- 
tons with use of the estimators in Eqs. H5.16|l and H5.17|l . respectively. Moreover, 
using the estimator in Eq. I|5.12(l we reconstruct the matrix elements 

Cn,m = a{m\b{m\-^){^\n)a\n)b, (5.33) 

which reveal the coherence of the twin-beam state. Theoretically one should 
have 

c„,™ = (i-icnrr"- (5.34) 

The estimators have been numerically evaluated by applying the Gauss method 
for calculating the integral in Eq. (|5.12() , which results in a fast and sufficiently 
precise algorithm with use of just 150 evaluation points. 
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N N 

Figure 5.2: Probability distribution for the total number of photons of the twin 
beams in Eq. H5.21|l for average number of photons per beam n — 2 obtained 
using the estimator in Eq. H5.17|l . The oscillation of the total photon- number 
probability due to the perfect correlation of the twin beams has been recon- 
structed by simulating 10^ data samples with quantum efficiency = 0.9 (on 
the left), and 2 x 10^ data samples rj = 0.8 (on the right). The theoretical 
probability (thick solid line) is superimposed onto the result of the Monte-Carlo 
experiment; the latter is shown by the thin solid line. Notice the dramatic 
increase of errors (in gray shade) versus N and for smaller 77 (From Ref. ^17,). 

In Fig. I5.1l a Monte-Carlo simulation of the joint photon-number probabil- 
ity distribution is reported. The simulated values compare very well with the 
theoretical ones. In Ref. ^21 a careful analysis of the statistical errors has been 
done for various twin-beam states by constructing histograms of deviations of 
the results from different simulated experiments from the theoretical ones. In 
comparison to the customary two-LO tomography of Ref. |87| , where for t] — 1 
the statistical errors saturate for increasingly large n and m, here we have sta- 
tistical errors that are slowly increasing versus n and m. This is due to the fact 
that the range of the estimators in Eq. H5.1(i|l increases versus n and m. Overall 
we find that for any given quantum efficiency the statistical errors are generally 
slightly larger than those obtained with the two-LO method. The convenience 
of using a single LO then comes with its own price tag. 

By using the estimator in Eq. H5.17|l the probability distribution for the total 
number of photons TV of the twin beams has been also constructed (Fig. I5.2|l . 
Notice the dramatic increase of error bars versus N and for smaller rj. 

Finally, in Fig. 15.81 we report the results of the tomographic measurement 
of Cn,m defined in Eq. (|5.33|) . Because the reconstructed Cn,m is close to the 
theoretically expected value in Eq. H5.34|l . these reveal the purity of the twin 
beams, which cannot be inferred from the thermal diagonal distribution of Fig. 
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Figure 5.3: Tomographic reconstruction of the matrix elements Cn.m = 
a{m\b{Tn\'i>){'i/\n)a\n)b of the twin beams in Eq. H5.21|l for average number of 
photons per beam n = 2, obtained using the estimator in Eq. (|5.12(l . On the 
left we used 10^ simulated data samples and quantum efhciency r] = 0.9; on the 
right 3 X 10^ data samples and 77 = 0.8. The coherence of the twin-beam state 
is easily recognized as Cn,m varies little for n + rn ^ constant in Eq. H5.21|l 
has been chosen real] . For a typical comparison between theoretical and exper- 
imental matrix elements and their relative statistical errors, see results in Fig. 
0(From Ref. [HI). 

The first experimental results of a measurement of the joint photon-number 
probability distribution for a two-mode quantum state created by a nondegen- 
erate optical parametric amplifier has been presented in Ref. [SHj. In this 
experiment, however, the twin beams are detected separately by two balanced- 
homodyne detectors. A schematic of the experimental setup is reported in Fig. 
15.41 and some experimental results are reported in Fig. 15.51 As expected for 
parametric fluorescence, the experiment has shown a measured joint photon- 
number probability distribution that exhibited up to 1.9 dB of quantum corre- 
lation between the two modes, with thermal marginal distributions. 
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Figure 5.4: A schematic of the experimental setup. NOPA: non-degenerate op- 
tical parametric amplifier; LOs: local oscillators; PBS: polarizing beam splitter; 
LPFs: low-pass filters; BPF: band-pass filter; G: electronic amplifier. Electron- 
ics in the two channels are identical (From Ref. (86jV 




Figure 5.5: Left: Measured joint photon-number probability distribution for the 
twin-beam state with average number of photons per beam n = 1.5 and 400000 
samples. Right: marginal distribution for the signal beam for the same data. 
The theoretical distribution is also shown. Very similar results are obtained for 
the idler beam (From Ref. |86J). 
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Chapter 6 

Applications to quantum 
measurements 

In this chapter we review a number of apphcations of quantum tomography 
related to some fundamental tests in quantum mechanics. 

First, we report the proposal of Ref. PU] for testing the nonclassicality of 
quantum states by means of an operational criterion based on a set of quantities 
that can be measured experimentally with some given level of confidence, even 
in the presence of loss, noise, and less-than-unity quantum efficiency. 

Second, we report the experiment proposed in Ref. for testing quantum 
state reduction. The state-reduction rule is tested using optical homodyne to- 
mography by directly measuring the fidelity between the theoretically-expected 
reduced state and the experimental state. 

Finally, we review some experimental results obtained at the Quantum Op- 
tics Lab of the University of Naples [32] about the reconstruction of coherent 
signals, together with application to the estimation of the losses introduced by 
simple optical components. 

6.1 Measuring the nonclassicality of a quantum 
state 

The concept of nonclassical states of light has drawn much attention in quantum 
optics imiHHllHniEniinillMlESlElinSllMl- The customary definition of 
nonclassicality is given in terms of the P-function presented in Sec. 12.11 a 
nonclassical state does not admit a regular positive P-function representation, 
namely, it cannot be written as a statistical mixture of coherent states. Such 
states produce effects that have no classical analogue. These kinds of states are 
of fundamental relevance not only for the demonstration of the inadequacy of 
classical description, but also for applications, e.g., in the realms of information 
transmission and interferometric measurements (91i .92,, .95^ . 
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We are interested in testing the nonclassicality of a quantum state by means 
of a set of quantities that can be measured experimentaUy with some given level 
of confidence, even in the presence of loss, noise, and less-than-unity quantum 
efficiency. The positivity of the P-function itself cannot be adopted as a test, 
since there is no viable method to measure it. As proved in Sec. 14. II only the 
generalized Wigner functions of order s < 1 — r/^^ can be measured, 77 being 
the quantum efhciency of homodyne detection. Hence, through this technique, 
all functions from s = 1 to s = cannot be recovered, i.e., we cannot obtain 
the P-function and all its smoothed convolutions up to the customary Wigner 
function. For the same reason, the nonclassicality parameter proposed by Lee 
|41| . namely, the maximum s-parameter that provides a positive distribution, 
cannot be experimentally measured. 

Among the many manifestations of nonclassical effects, one finds squeezing, 
antibunching, even-odd oscillations in the photon-number probability, and neg- 
ativity of the Wigner function jEHl EOl EH ESI EH EHl Ei lEDl- Any of these 
features alone, however, does not represent the univocal criterion we are look- 
ing for. Neither squeezing nor antibunching provides a necessary condition for 
nonclassicality The negativity of the Wigner function, which is well ex- 

hibited by the Fock states and the Schrodinger-cat-like states, is absent for the 
squeezed states. As for the oscillations in the photon-number probability, some 
even-odd oscillations can be simply obtained by using a statistical mixture of 
coherent states. 

Many authors E3 El EH have adopted the non-positivity of the phase- 
averaged P-function F{I) = ^ J^^ d4>P{I^^'^e^'^) as the definition for a non- 
classical state, since F{I) < invalidates Mandel's semiclassical formula 
of photon counting, i.e., it does not allow a classical description in terms of a 
stochastic intensity. Of course, some states can exhibit a "weak" nonclassical- 
ity inn], namely, a positive F{I), but with a non-positive P-function (a relevant 
example being a coherent state undergoing Kerr- type self-phase modulation). 
However, from the point of view of the detection theory, such "weak" non- 
classical states still admit a classical description in terms of positive intensity 
probability F{I) > 0. For this reason, we adopt non-positivity of F{I) as the 
definition of nonclassicality. 

6.1.1 Single-mode nonclassicality 

The authors of Rcfs. 93, 94, 96 have pointed out some relations between F{I) 
and generalized moments of the photon distribution, which, in turn, can be used 
to test the nonclassicality. The problem is reduced to an infinite set of inequali- 
ties that provide both necessary and sufficient conditions for nonclassicality |94| . 
In terms of the photon-number probability p(n) — {n\p\n) of the state with den- 
sity matrix p, the simplest sufficient condition involves the following three-point 
relation ^ EEl 

B{n) = {n + 2)p{n)p{n + 2) - {n + l)[p{n + l)f < . (6.1) 
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Higher-order sufficient conditions involve five-, seven-, . . . , {2k + l)-point re- 
lations, always for adjacent values of n. It is sufficient that just one of these 
inequalities is satisfied in order to assure the negativity of F{I). Notice that for 
a coherent state B{n) = identically for all n. 

In the following we show that quantum tomography can be used as a powerful 
tool for performing the nonclassicality test in Eq. 1)6. Ifl . For less-than-unity 
quantum efficiency (rj < 1), we rely on the concept of a "noisy state" p^, wherein 
the effect of quantum efficiency is ascribed to the quantum state itself rather 
than to the detector. In this model, the effect of quantum efficiency is treated 
in a Schrodinger-like picture, with the state evolving from p to pr/, and with j] 
playing the role of a time parameter. Such lossy evolution is described by the 
master equation [ST] 

p 

dtp{t) ^ - [2ap(t)a^ - a^ap(t) - p(t)a^a] , (6.2) 

wherein p{t) = p^ with t — — Inry/F. 

For the nonclassicality test, reconstruction in terms of the noisy state has 
many advantages. In fact, for non-unit quantum efficiency 77 < 1 the tomo- 
graphic method introduces errors for p{n) which are increasingly large versus n, 
with the additional limitation that quantum efficiency must be greater than the 
minimum value 77 = 0.5. On the other hand, the reconstruction of the noisy- 
state probabilities Pr]{n) — {n\p^\n) does not suffer such limitations, and even 
though all quantum features are certainly diminished in the noisy-state descrip- 
tion, nevertheless the effect of non-unity quantum efficiency does not change 
the sign of the P-function, but only rescales it as follows: 

P{z) ^ P,{z) = l^P{z/v'/') . (6.3) 

Hence, the inequality (|6.1|) still represents a sufficient condition for nonclassi- 
cality when the probabilities p{n) = {n\p\n) are replaced with Pjj{n) — (n|p^|n), 
the latter being given by a Bernoulli convolution, as shown in Eq. (|2.22|) . When 
referred to the noisy-state probabilities (n) , the inequality in Eq. (|6.1|l keeps 
its form and simply rewrites as follows 

B„{n) = {n + 2)p„(n)p,(n + 2) - {n + l)[p^{n + 1)]^ < . (6.4) 

The quantities B{n) and Br,{n) are nonlinear in the density matrix. Then, 
they cannot be measured by averaging a suitable estimator over the homodyne 
data. Hence, in the evaluation of B(n) one has to reconstruct the photon- 
number probabilities p(n), using the estimator 7?.^[|n)(?i|](a;, (^) in Eq. (|3.44|) . 
The noisy-state probabilities Prj{n) are obtained by using the same estimator 
for rj — 1, namely without recovering the convolution effect of non-unit quan- 
tum efficiency. Notice that the estimator does not depend on the phase of the 
quadrature. Hence, the knowledge of the phase of the local oscillator in the ho- 
modyne detector is not needed for the tomographic reconstruction, and it can 
be left fluctuating in a real experiment. 
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Regarding the estimation of statistical errors, they are generally obtained by 
dividing the set of homodyne data into blocks, as shown in Sec. 13.3.11 However, 
in the present case, the nonlinear dependence on the photon number probability 
introduces a systematic error that is vanishingly small for increasingly larger sets 
of data. Therefore, the estimated value of B{n) is obtained from the full set of 
data, instead of averaging the mean value of the different statistical blocks. 



0^ 
CP 



Figure 6.1: Tomographic measurement of B{n) (dashed trace) with the respec- 
tive error bars (superimposed in gray-shade) along with the theoretical values 
(solid trace) for a Schrodinger cat state with average photon number n = 5 
(left) ; for a phase-squeezed state with n = 5 and fisq = sinh^ r — 3 squeezing 
photons (right). In both cases the quantum efhciency is 77 = 0.8 and the number 
of simulated experimental data is 10^ (From Ref. |3U| ) . 

In Figs. l6.lll6.2l some numerical results from Ref. PU] are reported, which are 
obtained by a Monte-Carlo simulation of a quantum tomography experiment. 
The nonclassicality criterion is tested either on a Schrodinger-cat state \ip{a)) oc 
(|a) -I- I — a)) or on a squeezed state \a,r) = D{a)S{r)\0), wherein |a), D{a), 
and S{r) denote a coherent state with amplitude a, the displacement operator 
D{a) — e""*"" and the squeezing operator S{r) = ^<'{o}'^-a^)l2 ^ respectively. 
Fig. 16.11 shows tomographically-obtained values of B{n), with the respective 
error bars superimposed, along with the theoretical values for a Schrodinger-cat 
state and for a phase-squeezed state (r > 0). For the same set of states the 
results for Brj{n) obtained by tomographic reconstruction of the noisy state are 
reported in Fig. 16.21 Let us compare the statistical errors that affect the of B{n) 
and Bri{n) on the original and the noisy states, respectively. In the first case 
the error increases with n, whereas in the second it remains nearly constant, 
albeit with less marked oscillations in B,i{n) than those in B{n). 

The nonclassicality of the states here analyzed is experimentally verifiable, 
as Bjj{0) < by more than five standard deviations. In contrast, for coherent 
states one obtains small statistical fluctuations around zero for all n. Finally, 
we remark that the simpler test of checking for antibunching or oscillations in 
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Figure 6.2: Same as Fig. 16.11 but here for B,^{n) (From Ref. PU)'). 

the photon-number probability in the case of the phase-squeezed state (left of 
Figs. 16. ll and 16.2(1 would not reveal the nonclassical features of such a state. 

6.1.2 Two- mode nonclassicality 

In Ref. (30| it is also shown how quantum homodyne tomography can also 
be employed to test the nonclassicality of two-mode states. For a two-mode 
state nonclassicality is defined in terms of non-positivity of the following phase- 
averaged two-mode P-function j96| : 

F{h,l2,cP) = ^^'"d0iP(/y'e'*S/2'^V(^^+*)) . (6.5) 

In Ref. it is also proved that a sufficient condition for nonclassicality is 

C = ((ni - n2)2) - ((m - n2»2 - (m + na) < , (6.6) 

where ni and n2 are the photon-number operators of the two modes. 

A tomographic test of the inequality in Eq. (|6.6|) can be performed by aver- 
aging the estimators for the involved operators using Table ElTl Again, the value 
r] — 1 can be used to reconstruct the ensemble averages of the noisy state p^. 
As an example, we consider the twin-beam state of Eq. (|5.21|l . The theoretical 
value of C is given by C -2|^|V(1 - \^\^) < 0. With regard to the effect of 
quantum efficiency 77 < 1, the same argument still holds as for the single-mode 
case: one can evaluate C,, for the twin beams degraded by the effect of loss, and 
use 77 = 1 in the estimators. In this case, the theoretical value of Crj is simply 
rescaled, namely 

C,^~2rf\e/{l~m- (6.7) 

In Fig. 16.31 we report C,, vs. 1 — 77, with 77 ranging from 1 to 0.3 in steps 
of 0.05, for the twin beam in Eq. ((5.21() with = 0.5, corresponding to a 
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total average photon number (ni +712) = 2. The values of Cr/ result from a 
Monte-Carlo simulation of a homodyne tomography experiment with a sample 
of 4 X 10^ data. The nonclassicality test in terms of the noisy state gives values 
of Cr, that are increasingly near the classically positive region for decreasing 
quantum efficiency 77. However, the statistical error remains constant and is 
sufficiently small to allow recognition of the nonclassicality of the twin beams 
up to ?7 = 0.3. 
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Figure 6.3: Tomographic measurement of the nonclassical parameter C,; for twin 
beams in Eq. 15.2111 with = 0.5. The results are shown for different values 
of the quantum efficiency 77 (in steps of 0.05), and for each value the number of 
simulated data is 4 x lO'^. Statistical errors are shown in the gray shade (From 
Ref. I^). 

We conclude that quantum homodyne tomography allows one to perform 
nonclassicality tests for single- and two-mode radiation states, even when the 
quantum efhciency of homodyne detection is rather low. The method involves 
reconstruction of the photon- number probability or of some suitable function of 
the number operators pertaining to the noisy state, namely, the state degraded 
by the less-than-unity quantum efficiency. The noisy-state reconstruction is 
affected by the statistical errors; however, they are sufficiently small that the 
nonclassicality of the state can be tested even for low values of 77. For the cases 
considered here, we have shown that the nonclassicality of the states can be 
proved (deviation from classicality by many error bars) with lO'^-lO^ homodyne 
data. Moreover, since the knowledge of the phase of the local oscillator in the 
homodyne detector is not needed for the tomographic reconstruction, it can be 
left fluctuating in a real experiment. 

6.2 Test of state reduction 

In quantum mechanics the state reduction (SR) is still a very discussed rule. 
The so-called "projection postulate" was introduced by von Neumann j2] to ex- 
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plain the results from the Compton-Simons experiment, and it was generalized 
by Liiders |101j for measurements of observables with degenerate spectrum. The 
consistency of the derivation of the SR rule and its validity for generic measure- 
ments have been analyzed with some criticism |l(J2j . In a very general context, 
the SR rule was derived in a physically consistent way from the Schrodinger 
equation for the composite system of object and measuring apparatus |1U3| . An 
experiment for testing quantum SR is therefore a very interesting matter. Such 
a test in general is not equivalent to a test of the repeatability hypothesis since 
the latter holds only for measurements of observables that are described by 
self-adjoint operators. For example, joint measurements like the Arthurs-Kelly 
are not repeatable, as the reduced states are coherent states, which are not 
orthogonal. 

Quantum optics offers a possibihty of testing the SR, because several ob- 
servables can be chosen to perform different measurements on a fixed system. 
For instance, one can decide to perform either homodync or heterodyne, or 
photon-number detection. This is a unique opportunity; in contrast, in parti- 
cle physics the measurements are mostly quasi-classical and restricted to only 
a few observables. In addition, optical homodyne tomography allows a precise 
determination of the quantum system after the SR. 

A scheme for testing the SR could be based on tomographic measurements 
of the radiation density matrix after nondemolition measurements. However, 
such a scheme would reduce the number of observables that are available for 
the test. Instead, one can take advantage of the correlations between the twin 
beams of Eq. H5.21|l produced by a non-degenerate optical parametric amplifier 
(NOPA), in which case one can test the SR even for demolitive-type measure- 
ments. Indeed, if a measurement is performed on one of the twin beams, the SR 
can be tested by homodyne tomography on the other beam. This is precisely the 
scheme for an experimental test of SR proposed in Ref. |31| . which is reviewed 
in the following. 

The scheme for the SR test is given in Fig. 16.41 Different kinds of measure- 
ments can be performed on beam 1, even though here the SR only for heterodyne 
detection and photon-number detection will be considered. 

For a system described by a density operator p, the probability p{X)dX that 
the outcome of a quantum measurement of an observable is in the interval 
[A, A + dX) is given by Born's rule p{X)dX = Tr[p IlxdX] , where Ha is the POVM 
pertaining to the measurement that satisfies Ha > and / dX Il\ ~ I. For 
an exact measurement of an observable, which is described by a self- adjoint 
operator, IIa is just the projector over the eigenvector corresponding to the 
outcome A. In the case of the photon number a' a the spectrum is discrete and 
the POVM is ![„ = \m){m\ for integer eigenvalue m. For the Arthurs-Kelly 
joint measurement of the position and momentum (corresponding to a joint 
measurement of two conjugated quadratures of the field) we have the coherent- 
state POVM n„ = 7r-i|a)(a|. 

When on beam 1 we perform a measurement described by IIa, the reduced 
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Figure 6.4: Schematic of the proposed scheme for testing the SR for heterodyne 
detection. A NOPA generates a pair of twin beams (1 and 2). After heterodyn- 
ing beam 1, the reduced state of beam 2 is analyzed by homodyne tomography, 
which is conditioned by the heterodyne outcome. In place of the heterodyne 
detector one can put any other kind of detector for testing the SR on different 
observables. We also consider the case of direct photodetection (From Ref.|31p. 



normalized state of beam 2 is 

. Tri[|O(g|(n;,0i)] _ sni;st 

'^^^ - Tri,.[ie)(ei(n.«i)] - -wr ' ^'-'^ 

where O"^ denotes the transposed operator (on a fixed basis), S = (1— ICP)^^^'?'^ 
andp(A) = Tri_2[2n5|S^^] is the probability density of the measurement outcome 
A. In the limit of infinite gain |^| 1 one has p(A) oc IT^. For example, for 
heterodyne detection with outcome a, we have p{a) — \a*){a*\. 

If the readout detector on beam 1 has quantum efficiency rjr, Eq. ()5.8|l is 
replaced with 

where p^''{X) = Tri,2[S(n^'')'^S^], and n^"^ is the POVM for measurement with 
quantum efficiency rjr- As shown in Sec. 12.41 for heterodyne detection one has 
the Gaussian convolution 

'2. 



1 f d^Z 



with A^^ = (1 — rjr) /rjr. For direct photodetection Ilm = \m){m\ is replaced 
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with the Bernoulli convolution 



oo 



^"^-JliL) '^"(i - vry~"'\j)(j\ ■ (6-11) 

j—m 

The experimental test proposed here consists of performing conditional ho- 
modyne tomography on beam 2, given the outcome A of the measurement on 
beam 1. We can directly measure the "fidelity of the test" 

F{X)=Ty[p^^{X)p,^,U^)], (6.12) 

where p^''{X) is the theoretical state in Eq. (|6.9|l . and Pmcas(A) is the experimen- 
tally measured state on beam 2. Notice that we use the term "fidelity" even 
if F{X) is a proper fidelity when at least one of the two states is pure, which 
occurs in the limit of unit quantum efficiency r]r. In the following we evaluate 
the theoretical value of F{X) and compare it with the tomographic measured 
value. 

The fidelity (|6.12() can be directly measured by homodyne tomography using 
the estimator for the operator p^''{X), namely 



F{X)^ — / da;p,„(a;,<^;A)7^,Jp'"-(A)](x,^) , (6.13) 

where pri^{x, ip; A) is the conditional homodyne probability distribution for out- 
come A at the readout detector. 

For heterodyne detection on beam 1 with outcome a € C, the reduced state 
on beam 2 is given by the displaced thermal state 

p^^ia) = ViDi^){l - Vif^DHl) , (6.14) 

where 

^^^l + {r,r-m\\ 7=— «*• (6.15) 
The estimator in Eq. (|6.13|l is given by 

7^,Jp'M«)](-,^) = /^<i>^l,|;-;J^^(--7.)^) , (6.16) 
2rih - 775 V 2 277/1 - 775 J 

where 7,^ = Re(7e~"^), and ^{a,b;z) denotes the customary confluent hyper- 
geometric function. The estimator in Eq. H6.16|l is bounded for rjii > ^77^, then 
one needs to have 

Vh>\[l-\ea'Vr)] ■ (6.17) 

As one can see from Eq. (|6.17|l . for rjh > 0.5 the fidelity can be measured for any 
value of rjr and any gain parameter ^ of the NOPA. We recall that the condition 
•qh > 0.5 is required for the measurement of the density matrix. However, in 
this direct measurement of the fidelity, the reconstruction of the density matrix 
is bypassed, and we see from Eq. H6.17|l that the bound rjh ~ 0.5 can be lowered. 
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The measured fidelity F{a) in Eq. Ht).13|l with p^''{a) as given in Eq. H6.14|l 
must be compared with the theoretical value 



Fa,^Vi/{2-Vi) , (6.18) 

that is independent of a. 

For direct photodetection on beam 1 with outcome n, the reduced state on 
beam 2 is given by 

The estimator for the fidelity measurement is 



$ 1, _; ^^iliS I^a;^ (B.20) 

2r]h - ri^ + z \ 2 277^-77^+2 / 

We see that the same bound of Eq. (|6.17|) holds. In this case the measured 
fidelity i^(77) must be compared with the theoretical value 

= vl^^'^F (277 + 1, 277 + 1; 1; (1 - 77^)2) , (6.21) 

where -^(0, 6; c; z) denotes the customary hypergeometric function. 

Several simulations have been reported in Ref. [2] for both heterodyne 
and photodetection on beam 1. In the former case the quadrature probability 
distribution has been simulated, pertaining to the reduced state H6.14II on beam 
2, and averaged the estimators in Eq. (|6.16|) . In the latter case the reduced 
state (|6.19|) and the estimators in Eq. (|6.20(l have been used. 

Numerical results for the fidelity was thus obtained for different values of 
the quantum efficiencies 77^ and 77;^, and of the NOPA gain parameter ^. A deci- 
sive test can be performed with samples of just a few thousand measurements. 
The statistical error in the measurement was found rather insensitive to both 
quantum efficiencies and NOPA gain. 



6.3 Tomography of coherent signals and appU- 
cations 

Quantum homodyne tomography has been proved useful in various experimental 
situations, such as for measuring the photon statistics of a semiconductor laser 
for determining the density matrix of a squeezed vacuum [TT] and the 
joint photon-number probability distribution of a twin beam created by a non- 
degenerate optical parametric amplifier [5^1 , and, finally, for reconstructing the 
quantum states of spatial modes with an array detector |104| . In this section we 
review some experimental results about homodyne tomography with coherent 
states, with application to the estimation of the loss introduced by simple optical 
components [?2] . 
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The experiment has been performed in the Quantum Optics Lab of the Uni- 
versity of Naples, and its the schematic is presented in Fig. 16.51 The principal 
radiation source is provided by a monolithic Nd:YAG laser (w 50 mW at 1064 
nm; Lightwave, model 142). The laser has a linewidth of less than 10 kHz/ms 
with a frequency jitter of less than 300 kHz/s, while its intensity spectrum is 
shot-noise limited above 2.5 MHz. 



FR HIVP, HWP, 




Figure 6.5: Schematic of the experimental set-up. A Nd:YAG laser beam is 
divided into two beams, the first acting as a strong local oscillator, the second 
representing the signal beam. The signal is modulated at frequency with a 
defined modulation depth to control the average photon number in the generated 
coherent state. The tomographic data are collected by a homodyne detector 
whose difference photocurrent is demodulated and then acquired by a digital 
oscilloscope (From Ref. |32)'). 

The laser emits a linearly polarized beam in a TEMOO mode, which is split 
in two parts by a beam splitter. One part provides the strong local oscillator 
for the homodyne detector. The other part, typically less than 200 /iW, is the 
homodyne signal. The optical paths traveled by the local oscillator and the 
signal beams are carefully adjusted to obtain a visibility typically above 75% 
measured at one of the homodyne output port. The signal beam is modulated, 
by means of a phase electro-optic modulator (EOM, Linos Photonics PM0202), 
at 4MHz, and a halfwave plate (HWP2, HWP3) is mounted in each path to 
carefully match the polarization state at the homodyne input. 

The detector is composed by a 50-;-50 beam splitter (BS), two amplified 
photodiodes (PDl, PD2), and a power combiner. The difference photocurrent 
is demodulated at 4MHz by means of an electrical mixer. In this way the 
detection occurs outside any technical noise and, more important, in a spectral 
region where the laser does not carry excess noise. 

The phase modulation added to the signal beam moves a certain number of 
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photons, proportional to the square of the modulation depth, from the carrier 
optical frequency uj to the side bands at (jj ± f2 so generating two weak coherent 
states with engineered average photon number at frequencies w ± fJ. The sum 
sideband modes is then detected as a controlled perturbation attached to the 
signal beam. The demodulated current is acquired by a digital oscilloscope 
(Tektronix TDS 520D) with 8 bit resolution and record length of 250000 points 
per run. The acquisition is triggered by a triangular shaped waveform applied 
to the PZT mounted on the local oscillator path. The piezo ramp is adjusted 
to obtain a 2tt phase variation between the local oscillator and the signal beam 
in an acquisition window. 

The homodyne data to be used for tomographic reconstruction of the state 
have been calibrated according to the noise of the vacuum state. This is obtained 
by acquiring a set of data leaving the signal beam undisturbed while scanning 
the local oscillator phase. It is important to note that in case of the vacuum 
state no role is played by the visibility at the homodyne beam-splitter. 

The tomographic samples consist of N homodyne data {xj, ^j}j=i,...,N with 
phases ipj equally spaced with respect to the local oscillator. Since the piezo 
ramp is active during the whole acquisition time, we have a single value Xj for 
any phase ipj. From calibrated data we first reconstruct the quantum state 
of the homodyne signal. According to the experimental setup, we expect a 
coherent signal with nominal amplitude that can be adjusted by varying the 
modulation depth of the optical mixer. However, since we do not compensate 
for the quantum efficiency of photodiodes in the homodyne detector (77 ~ 90%) 
we expect to reveal coherent signals with reduced amplitude. In addition, the 
amplitude is further reduced by the non-maximum visibility (ranging from 75% 
to 85%) at the homodyne beam-splitter. 

In Fig. It). 51 we report a typical reconstruction, together with the reconstruc- 
tion of the vacuum state used for calibration. For both states, we report the raw 
data, the photon number distribution /)„„, and a contour plot of the Wigner 
function. The matrix elements are obtained by sampling the corresponding es- 
timators in Eq. H3.44|l . whereas the confidence intervals for diagonal elements 
are given by 5pnn = Ap/ a/^V, Ap being the rms deviation of the estimator over 
data. For off-diagonal elements the confidence intervals are evaluated for the 
real and imaginary part separately. 

In order to see the quantum state as a whole, we also report the reconstruc- 
tion of the Wigner function of the field, which can be expressed in terms of the 
matrix elements as the discrete Fourier transform 

00 00 

I^(a,a*) = Re^e^'^'^^A(n,d;|a|)p„,„+d (6.22) 

d=0 n=a 

where ip — arga, and 



A(n,d;|«|) = (-)"2(2-M|2«^/7;^e-2|"I^L:^(|2an , (6.23) 
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L^(a;) denoting the Laguerre polynomials. Of course, the series in Eq. H6.22|l 
should be truncated at some point, and therefore the Wigner function can be 
reconstructed only at some finite resolution. 
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Figure 6.6: Reconstruction of the quantum state of the signal, and of the vacuum 
state used for calibration. For both states, from left to right, we report the raw 
data, a histogram of the photon number distribution, and a contour plot of 
the Wigner function. The reconstruction has been performed by a sample of 
N = 242250 homodyne data. The coherent signal has an estimated average 
photon number equal to (a^a) = 8.4. The solid line denotes the theoretical 
photon distribution of a coherent state with such number of photons. Statistical 
errors on matrix elements are about 2% . The slight phase asymmetry in the 
Wigner distribution corresponds to a value of about 2% of the maximum (From 
Ref. El). 

Once the coherence of the signal has been established we may use homodyne 
tomography to estimate the loss imposed by a passive optical component like 
an optical filter. The procedure may be outlined as follows. Wc first estimate 
the initial mean photon number uq — |aoP of the signal beam, and then the 
same quantity inserting an optical neutral density filter in the signal path. If F 
is the loss parameter, then the coherent amplitude is reduced to ar = aoe~^ , 
and the intensity to fir = nQe~'^^ . 

The estimation of the mean photon number can be performed adaptively on 
data, using the general method presented in Sec. 13.4.21 One takes the average 
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of the estimator 

n[a)a] {x, ip) = 2x^ -]^+ fie'^^ + fi*e-'^^ , (6.24) 

where /x is a parameter to be determined in order to minimize fluctuations. As 
proved in Ref. [22 one has fj, = —l/2{a^'^), which itself can be obtained from 
homodyne data. In practice, one uses the data sample twice: first to evaluate 
/X, then to obtain the estimate for the mean photon number. 
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Figure 6.7: Estimation of the mean photon number of a coherent signal as a 
function of the loss imposed by an optical filter. Three set of experiments, 
corresponding to three different initial amplitudes are reported. Open circles 
are the tomographic determinations, whereas solid line denotes the expected 
values, as follow from nominal values of loss and visibility at homodyne detector. 
Statistical errors are within the circles (From Ref. |32p. 

In Fig. 16.71 the tomographic determinations of nr are compared with the 
expected values for three set of experiments, corresponding to three different 
initial amplitudes. The expected values are given by nr = nQe~^^V, where 
r is the value obtained by comparing the signal dc currents Iq and Ir at the 
homodyne photodiodes and V = Vr/Vo is the relative visibility. The solid line 
in Fig. 16.71 denotes these values. The line is not continuous due to variations 
of visibility. It is apparent from the plot that the estimation is reliable in the 
whole range of values we could explore. It is worth noting that the estimation is 
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absolute, i.e. it does not depend on the knowledge of the initial amplitude, and 
it is robust, since it can be performed independently on the quantum efficiency 
of the homodyne detector. 

One may notice that the estimation of loss can be pursued also by measur- 
ing an appropriate observable, typically the intensity of the light beam with 
and without the filter. However, this is a concrete possibility only for high 
amplitude signals, whereas losses on weak coherent states cannot be properly 
characterized neither by direct photocounting using photodiodes (due to the 
low quantum efficiency and large fluctuations), nor by avalanche photodetec- 
tors (due to the impossibility of discriminating among the number of photons). 
On the contrary, homodyne tomography provides the mean intensity (actually 
the whole photon distribution) independently on the signal level, thus allowing 
a precise characterization also in the quantum regime. Indeed, in Ref. |22| 
adaptive tomographic determination of the mean photon number has been ex- 
tensively applied to (numerically simulated) homodyne data for coherent states 
of various amplitudes. The analysis has shown that the determination is reliable 
also for small samples and that precision is not much affected by the intensity 
of the signal. 
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Chapter 7 

Tomography of a quantum 
device 

If we want to determine experimentally the operation of a quantum device, we 
need, by definition, quantum tomography. In fact, the characterization of the 
device operation could be done by running a basis of possible known inputs, and 
determining the corresponding outputs by quantum tomography. In quantum 
mechanics the inputs are density operators, and the role of the transfer matrix 
is played by the so-called quantum operation of the device, here denoted by £. 
Thus the output state pout (a part from a possible normalization) is given by 
the quantum operation applied to the input state as follows 

Pout = £{pin)- (7.1) 

Since the set of states p actually belongs to a space of operators, this means that 
if we want to characterize £ completely, we need to run a complete orthogonal 
basis of quantum states \n) (n = 0, 1, 2, . . .), along with their linear combina- 
tions "^d"-') + i'^l'iT'")), with k — 0, 1,2,3 and i denoting the imaginary unit. 
However, the availability of such a set of states in the lab is, by itself, a very 
hard technological problem. For example, for an optical device, the states \n) 
are those with a precise number n of photons, and, a part from very small n — 
say at most n=2 — they have never been achieved in the lab, whereas preparing 
their superpositions remains a dream for experimentalists, especially if n ^ 1 
(a kind of Schrodinger kitten states). 

The idea of achieving the quantum operation of a device by scanning the 
inputs and making tomography of the corresponding output is the basis of the 
early methods proposed in Refs. jl()5l ll()()j . Due to the mentioned problems 
in the availability of input states, both methods have limited application. The 
method of Ref. |1U5| has been designed for NMR quantum processing, whereas 
the method of Ref. |106| was conceived for determining the Liouvillian of a 
phase-insensitive amplifier, namely for a case in which the quantum operation 
has no off-diagonal matrix elements, to evaluate which one needs the superpo- 
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sitions -^d^') + with k — 0, 1,2,3 mentioned above. The problem of 

availabihty of input spates and their superpositions was partiaUy solved by the 
method of Ref. |1U7) , where it was suggested to use randomly drawn coherent 
states to estimate the quantum operation of an optical device via a maximum 
likelihood approach. This method, however, cannot be used for quantum sys- 
tems different from the em radiation — such as finite dimensional systems, i.e. 
qubits — due to the peculiarity of coherent states. The solution to the problem 
came with the recent method of Ref. , where the problem of the availability 
of input states was solved by using a single bipartite entangled input, which 
is equivalent to run all possible input states in a kind of "quantum parallel" 
fashion (bipartite entangled states are nowadays easily available in most quan- 
tum systems of interest). The method is also very simple and effective, and 
its experimental feasibility (for single-photon polarization-encoded qubits) has 
been already demonstrated in a recent experiment performed in the Francesco 
De Martini laboratory in Roma La Sapienza jl(J8j . In the next sections we will 
review the general method and report some computer simulated results from 
Ref. 



7.1 The method 

As already mentioned, the most description of a general state-transformation in 
quantum mechanics is given in terms of the so-called quantum operation. The 
state transformation due to the quantum operation £ is given as follows 

" ■ (7.2) 



" Tr(f (p)) • 

The transformation occurs with probability given by p = Tt[£{p)] < 1. The 
quantum operation £ is a linear, trace-decreasing completely positive (CP) map. 
We remind that a map is completely positive if it preserves positivity generally 
when applied locally to an entangled state. In other words, upon denoting by 
X the identical map on the Hilbert space /C of a second quantum system, the 
extended map £ ®T onli.® K, is positive for any extension /C. Typically, the 
CP map is written using a Kraus decomposition |109| as follows 

£{p)=Y,KnpKi, (7.3) 

n 

where the operators satisfy 

Y.KiKa<I. (7.4) 

n 

The transformation 1)7. 3|l occurs with generally non-unit probability Tr[£(/9)] < 
1, and the probability is unit independently on p when £ is trace-preserving, 
i.e. when we have the equal sign in Eq. H7.4|l . The particular case of unitary 
transformations corresponds to having just one term Ki = U in the sum 1)7. 3|l . 
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with U unitary. However, one can consider also non- unitary operations with 
one term only, namely 

£{p) = ApA^ , (7.5) 

where A is a contraction, i. e. \\A\\ < 1. Such operations leave pure states as 
pure, and describes, for example, the state reduction from a measurement appa- 
ratus for a particular fixed outcome that occurs with probability Tr[pA^^yl] < 1. 

In the following we will use the notation for bipartite pure states introduced 
in Eq. (|2.45|) . and we will denote by O'^ and O* the transposed and the conjugate 
operator of O with respect to some pre-chosen orthonormal basis. 

The basic idea of the method in Ref. |2S] is the following. An unknown 
quantum operation £ can be determined experimentally through quantum to- 
mography, by exploiting the following one-to-one correspondence £ <-> Rg be- 
tween quantum operations £ and positive operators on two copies of the 
Hilbert space T-L^H. 

Rs ^£(^I{\I)){{I\), £{p)^Tt2[I®P^R£]. (7.6) 

Notice that the vector |/)) represents a (unnormalized) maximally entangled 
state. If we consider a bipartite input state lip)) and operate with £ only on one 
Hilbert space as in Fig. 17.11 the output state is given by 

R{^/j)=£^imU\). (7.7) 

For invertible tp the two matrices R{I) = Rs and R{ip) are related as follows 

R{I) = (/ ® V"^^^(V')(/ ® V-"^*) ■ (7.8) 

Hence, the (four-index) quantum operation matrix Rg can be obtained by esti- 
mating via quantum tomography the following ensemble averages 

{{i,j\R{I)\l,k))^TT[R{^l;){\l)m®i;-'*\k){M-'*)] . (7.9) 

Then one simply has to perform a quantum tomographic estimation, by mea- 
suring jointly two observables Xx and X'^^ from two quorums {Xa} and {^^^} 
for the two entangled quantum systems. 

7.2 An example in the optical domain 

In Ref. |l25j it is shown that the proposed method for quantum tomography 
of a device can be actually performed using joint homodyne tomography on a 
twin-beam from downconversion of vacuum, with an experimental setup similar 
to that used in the experiment in Ref. The feasibility analysis considers, 

as an example, the experimental determination of the quantum operation corre- 
sponding to the unitary displacement operator D{z) — e^"^"^ °. The pertaining 
matrix R{I) is given by 

RiI)^\Diz)))({Diz)\, (7.10) 
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Figure 7.1: General scheme of the method for the tomographic estimation of a 
quantum operation. Two identical quantum systems are prepared in a bipartite 
state with invertible ^. One of the two systems undergoes the quantum 
operation £, whereas the other is left untouched. At the output one performs 
a quantum tomographic estimation, by measuring jointly two observables X\ 
and X'^ from two quorums {^a} and {X'^} for the two Hilbert spaces, such 
as two different quadratures of the two field modes in a two-mode homodyne 
tomography (From Ref. |25|). 



which is the (unnormalizable) eigenstate of the operator a — with eigenvalue 
z, as shown in Sec. 12.41 As an input bipartite state, one uses the twin-beam 
from parametric downconversion of Eq. H5.21|l . which is clearly invertible, since 

V' = VT~^r'", V"' = ^^^=r"'" ■ (7.11) 
^ ' ' ^ ^1- ie|2) ^ ' 

The experimental apparatus is the same as in the experiment of Ref. |86) . 
where the twin-beam is provided by a non-degenerate optical parametric am- 
plifier (a KTP crystal) pumped by the second harmonic of a Q-switched mode- 
locked Nd:YAG laser, which produces a 100-MHz train of 120-ps duration pulses 
at 1064 nm. The orthogonally polarized twin beams emitted by the KTP crystal 
(one of which is displaced of D{z) by a nearly transparent beam splitter with a 
strong local oscillator) are separately detected by two balanced homodyne de- 
tectors that use two independent local oscillators derived from the same laser. 
This provides the joint tomography of quadratures X^i (g) X^n needed for the 
reconstruction. The only experimental problem which still need to be addressed 
(even though is practically solvable) with respect to the original experiment of 
Ref. is the control of the quadrature phases 0' and 0" with respect to the 
LO, which in the original experiment were random. 

In Fig. 17.21 the results of a simulated experiment are reported, for displace- 
ment parameter z — 1, and for some typical values of the quantum efficiency t] 
at homodyne detectors and of the total average photon number n of the twin 
beam. The diagonal elements Ann ~ {n\D{z)\n) — [(ri|(n|i?£)(2)|n)n)]^/^ are 
plotted for the displacement operator with z = 1 . The reconstructed values are 
shown by thin solid line on an extended abscissa range, with their respective 
error bars in gray shade, and compared to the theoretical probability (thick 
solid line). A good reconstruction of the matrix can be achieved in the given 
range with n 1, quantum efficiency as low as 77 = 0.7, and 10^ 10^ data. 
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The number of data can be decreased by a factor 100 — 1000 using the tomo- 
graphic max-hkehhood techniques of Ref. however at the expense of the 
complexity of the algorithm. Improving quantum efficiency and increasing the 
amplifier gain (toward a maximally entangled state) have the effect of making 
statistical errors smaller and more uniform versus the photon labels n and m of 
the matrix A„m- 



in 
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Figure 7.2: Homodyne tomography of the quantum operation corresponding 
to the unitary displacement operator D[z), with z — \. The reconstructed 
diagonal elements Ann — are shown (thin solid line on an extended 

abscissa range, with their respective error bars in gray shade) , compared to the 
theoretical value (thick solid line). Similar results are obtained for off-diagonal 
terms. The reconstruction has been achieved using at the input the twin beam 
state of Eq. (|5.21|l . with total average photon number n and quantum efficiency 
at homodyne detectors 77. Left: n = 5, t] = 0.9, and 150 blocks of 10^ data have 
been used. Right: n — 3, rj — 0.7, and 300 blocks of 2 • 10^ data have been used 
(From Ref. 



It is worth emphasizing that the quantum tomographic method of Ref. [53] 
for measuring the matrix of a quantum operation can be much improved by 
means of a max-likelihood strategy aimed at the estimation of some unknown 
parameters of the quantum operation. In this case, instead of obtaining the ma- 
trix elements of R{I) from the ensemble averages in Eq. (|7.9|l . one parametrizes 
R{I) in terms of unknown quantities to be experimentally determined, and the 
likelihood is maximized for the set of experimental data at various randomly 
selected (tensor) quorum elements, keeping the same fixed bipartite input state. 
This method is especially useful for a very precise experimental comparison be- 
tween the characteristics of a given device (e.g. the gain and loss of an active 
fiber) with those of a quantum standard reference. 
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Chapter 8 



Maximum-likelihood 
method in quantum 
estimation 



Quantum estimation of states, observablcs and parameters is, from very basic 
principles, matter of statistical inference from a population sampling, and the 
most comprehensive quantum estimation procedure is quantum tomography. As 
we have shown in Chapter 3, the expectation value of an operator is obtained by 
averaging an estimator over the experimental data of a "quorum" of observables. 
The method is very general and efhcient, however, in the averaging procedure, 
we have fluctuations which result in relatively large statistical errors. 

Another relevant strategy, the maximum-likelihood (ML) method, can be 
used for measuring unknown parameters of transformation on a given state 

or for measuring the matrix elements of the density operator itself 
The ML strategy [TTUl [TTT] is an entirely different approach to quantum state 
measurement compared to the standard quantum-tomographic techniques. The 
ML procedure consists in finding the quantum state, or the value of the pa- 
rameters, that are most likely to generate the observed data. This idea can be 
quantified and implemented using the concept of the likelihood functional. 

As regards state estimation, the ML method estimates the quantum state 
as a whole. Such a procedure incorporates a priori knowledge about relations 
between elements of the density matrix. This guarantees positivity and normal- 
ization of matrix, with the result of a substantial reduction of statistical errors. 
Regarding the estimation of specific parameters, we notice that in many cases 
the resulting estimators are efficient, unbiased and consistent, thus providing a 
statistically reliable determination. 

As we will show, by using the ML method only small samples of data are 
required for a precise determination. However, we want to emphasize that 
such method is not always the optimal solution of the tomographic problem, 
since it suffers from some major limitations. Besides being biased due to the 
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Hilbert space truncation — even though the bias can be very small if, from other 
methods, we know where to truncate — it cannot be generalized to the estimation 
of any ensemble average, but just of a set of parameters from which the density 
matrix depends. In addition, for increasing number of parameters the method 
has exponential complexity. 

In the following we will review the ML methods proposed in Refs. and 
|33|. by deriving the likelihood functional, and applying the ML method to 
the quantum state reconstruction, with examples for both radiation and spin 
systems, and, finally, considering the ML estimation for the relevant class of 
Gaussian states in quantum optics. 



8.1 Maximum likelihood principle 

Here we briefly review the theory of the maximum-likelihood (ML) estimation of 
a single parameter. The generalization to several parameters, as for example the 
elements of the density matrix, is straightforward. The only point that should 
be carefully analyzed is the parameterization of the multidimensional quantity 
to be estimated. In the next section the specific case of the density matrix will 
be discussed. 

Let p{x\X) the probability density of a random variable a:, conditioned to 
the value of the parameter A. The form of p is known, but the true value of 
A is unknown, and will be estimated from the result of a measurement of x. 
Let xi, a;2, xn be a random sample of size N. The joint probability density 
of the independent random variable xi,X2, ■■■,xn (the global probability of the 
sample) is given by 

C{xi,X2, ...,a;Ar|A) = nf^ip(a;fc|A) , (8.1) 

and is called the likelihood function of the given data sample (hereafter we will 
suppress the dependence of C on the data) . The maximum-likelihood estimator 
(MLE) of the parameter A is defined as the quantity Xmi = Xmi{{xk}) that 
maximizes C{X) for variations of A, namely Xmi is given by the solution of the 
equations 

dC{X) d^i^iX) 
The first equation is equivalent to dL/dX — where 

N 

L{X) - log£{X) = ^logp(xfelA) (8.3) 
fe=i 

is the so-called log-likelihood function. 

In order to obtain a measure for the confidence interval in the determination 
of Am; we consider the variance 

[Xr,ai{xk})~Xf . (8.4) 



]^c?Xfep(a;A;|A) 
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In terms of the Fisher information 



F 



it is easy to prove that 



dx 



dp{x\X) 
dX 



pix\X) 



(8.5) 



NF 



(8.6) 



where TV is the number of measurements. The inequahty in Eq. H8.6|l is known 
as the Cramer-Rao bound |112| on the precision of the ML estimation. Notice 
that this bound holds for any functional form of the probability distribution 
p{x\X), provided that the Fisher information exists VA and d\p{x\X) exists Vx. 
When an experiment has "good statistics" (i.e. for a large enough data sample) 
the Cramer-Rao bound is saturated. 



8.2 ML quantum state estimation 

In this section we review the method of the maximum likelihood estimation of 
the quantum state of Ref. [231 , focusing attention to the cases of homodyne and 
spin tomography. 

We consider an experiment consisting of N measurements performed on 
identically prepared copies of a given quantum system. Each measurement is 
described by a positive operator- valued measure (POVM). The outcome of the 
ith measurement corresponds to the realization of a specific element of the 
POVM used in the corresponding run, and we denote this element by Hi. The 
likelihood is here a functional of the density matrix C{p) and is given by the 
product 

N 

C{p) = l[TT{pIl^) , (8.7) 

i=l 

which represents the probability of the observed data. The unknown element 
of the above expression, which we want to infer from data, is the density ma- 
trix describing the measured ensemble. The estimation strategy of the ML 
technique is to maximize the likelihood functional over the set of the density 
matrices. Several properties of the likelihood functional are easily found, if we 
restrict ourselves to finite dimensional Hilbert spaces. In this case, it can be 
easily proved that £{p) is a concave function defined on a convex and closed 
set of density matrices. Therefore, its maximum is achieved either on a single 
isolated point, or on a convex subset of density matrices. In the latter case, the 
experimental data are insufficient to provide a unique estimate for the density 
matrix using the ML strategy. On the other hand, the existence of a single 
maximum allows us to assign unambiguously the ML estimate for the density 
matrix. 

The ML estimation of the quantum state, despite its elegant general formu- 
lation, results in a highly nontrivial constrained optimization problem, even if 
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we resort to purely numerical means. The main difBculty lies in the appropri- 
ate parameterization of the set of all density matrices. The parameter space 
should be of the minimum dimension in order to preserve the maximum of the 
likelihood function as a single isolated point. Additionally, the expression of 
quantum expectation values in terms of this parameterization should enable 
fast evaluation of the likelihood function, as this step is performed many times 
in the course of numerical maximization. 

For such purpose one introduces (23 a parameterization of the set of density 
matrices which provides an efficient algorithm for maximization of the likelihood 
function. We represent the density matrix in the form 

p = T^T , (8.8) 

which automatically guarantees that p is positive and Hermitian. The remaining 
condition of unit trace Trp = 1 will be taken into account using the method of 
Lagrange multipliers. In order to achieve the minimal parameterization, we 
assume that T is a complex lower triangular matrix, with real elements on the 
diagonal. This form of T is motivated by the Cholesky decomposition known 
in numerical analysis for arbitrary non negative Hermitian matrix. For 

an Af-dimensional Hilbert space, the number of real parameters in the matrix 
T is Af + 2M{M — l)/2 = NP, which equals the number of independent real 
parameters for a Hermitian matrix. This confirms that such parameterization 
is minimal, up to the unit trace condition. 

In numerical calculations, it is convenient to replace the likelihood functional 
by its natural logarithm, which of course does not change the location of the 
maximum. Thus the log-likelihood function subjected to numerical maximiza- 
tion is given by 

N 

L{T) = ^ In Tr(r^TH,) - ATr(T^r) , (8.9) 

i=l 

where A is a Lagrange multiplier accounting for normalization of p. Writing p in 
terms of its eigenvectors ji/ip) as p — J/^ I V'm ) (V'p 1 1 with real y^, the maximum 
likelihood condition dL/dy,y = reads 

N 

Ay. = ^[2/.(V.|H,|V^,)/Tr(pn,)] , (8.10) 

i=l 

which, after multiplication by i/i, and summation over v, yields X — N. The 
Lagrange multiplier then equals the total number of measurements N. 

This formulation of the maximization problem allows one to apply standard 
numerical procedures for searching the maximum over the real parameters of 
the matrix T. The examples presented below use the downhill simplex method 

The first example is the ML estimation of a single-mode radiation field. 
The experimental apparatus used in this technique is the homodyne detector. 
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According to Sec. 12.41 the homodync measurement is described by the positive 
operator-valued measure 



n{x;ip) 



2-q 



■ exp 



2r] 



1 



(8.11) 



7r(l-77) 

where rj is the detector efficiency, and )/2 is the quadrature 

operator at phase ip. 




Figure 8.1: Reconstruction of the density matrix of a single-mode radiation 
field by the ML method. The plot shows the matrix elements of a coherent 
state (left) with {a) a) — 1 photon, and for a squeezed vacuum (right) with 
(a^a) = 0.5 photon. A sample of 50000 simulated homodyne data for quantum 
efficiency rj = 80% has been used (From Ref. |23|). 

After N measurements, we obtain a set of pairs (x^; ipi), where i = 1, . . . , iV. 
The log-likelihood functional is given by Eq. (|8.9|l with 11^ = T-L{xi\Lpi). Of 
course, for a light mode it is necessary to truncate the Hilbert space to a fi- 
nite dimensional basis. We shall assume that the highest Fock state has M — \ 
photons, i.e. that the dimension of the truncated Hilbert space is M . For the ex- 
pectation Tr[TtrH(a;; Lp)] it is necessary to use an expression which is explicitly 
positive, in order to protect the algorithm against occurrence of small negative 
numerical arguments of the logarithm function. A simple derivation yields 

^{k\T\n + j)Bn+j,n{n\^x)i 



A/-1 k 

T,[T^Tn{x-p^)\ = ^Y.ll 



Tl = 
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Bn 



n + j 
n 



1/2 
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1/4 



\/2"n! 
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(8.14) 
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are the eigenstates of the harmonic oscillator in the position representation — 
Hn{x) being the nth Hermite polynomial. 

The ML technique can be applied to reconstruct the density matrix in the 
Fock basis from Monte Carlo simulated homodyne statistics. Fig. 18.11 depicts 
the matrix elements of the density operator as obtained for a coherent state 
and a squeezed vacuum, respectively. Remarkably, only 50000 homodyne data 
have been used for quantum efficiency rj = 80%. We recall that in quantum 
homodyne tomography the statistical errors are known to grow rapidly with 
decreasing efficiency rj of the detector (23 EDI- In contrast, the elements of the 
density matrix reconstructed using the ML approach remain bounded, as the 
whole matrix must satisfy positivity and normalization constraints. This results 
in much smaller statistical errors. As a comparison one could see that the same 
precision of the reconstructions in Fig. lS.ll could be achieved using 10^-10® data 
samples with conventional quantum tomography. On the other hand, in order to 
find numerically the ML estimate we need to set a priori the cut-off parameter 
for the photon number, and its value is limited by increasing computation time. 




Figure 8.2: ML reconstruction of the density matrix of a two-mode radiation 
field. On the left the matrix elements obtained for the state = (|00) + 
|11))/a/2; on the right for |*2) = (|01) + |10))/\/2. For we used 100000 
simulated homodyne data and rj = 80%; for |'I'2) we used 20000 data and 
T] = 90% (From Ref. 

Another relevant example is the reconstruction of the quantum state of two- 
mode field using single-LO homodyning of Chapter 5. Here, the full joint density 
matrix can be measured by scanning the quadratures of all possible linear com- 
binations of modes. For two modes the measured quadrature operator is given 

by 

X{9,%,iJi) = ^{ae-''l'° cose + be-''!'^ sine + h.c.) , (8.15) 

where {9, ipQjipi) £ S^x[0, 2tt\ , S*^ being the Poincare sphere and one phase rang- 
ing between and 2n. In each run these parameters are chosen randomly. The 
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POVM describing the measurement is given by the right-hand side of Eq. (|8.11|) , 
with X^p replaced by X{0, ipo,tpi). An experiment for the two orthogonal states 
l^fi) = (|00) + |11))/V2 and l^a) = (|01) + |10))/V2 has been simulated, in 
order to reconstruct the density matrix in the two-mode Fock basis using the 
ML technique. The results are reported in Fig. 18.21 

The ML procedure can also be applied for reconstructing the density matrix 
of spin systems. For example, let us consider TV repeated preparations of a 
pair of spin-1/2 particles. The particles are shared by two parties. In each run, 
the parties select randomly and independently from each other a direction along 
which they perform a spin measurement. The obtained result is described by the 
joint projection operator (spin coherent states \llb\ ) !Fi = jfi^, fjf )(ri^, fif |, 
where and Of are the vectors on the Bloch sphere corresponding to the 
outcomes of the ith run, and the indices A and B refer to the two particles. As 
in the previous examples, it is convenient to use an expression for the quantum 
expectation value Tr(T^^J^i) which is explicitly positive. The suitable form is 



where is an orthonormal basis in the Hilbert space of the two particles. The 
result of a simulated experiment with only 500 data for the reconstruction of 
the density matrix of the singlet state is shown in Fig. 18.31 



Figure 8.3: ML reconstruction of the density matrix of a pair of spin-1/2 parti- 
cles in the singlet state. The particles are shared by two parties. In each run, 
the parties select randomly and independently from each other a direction along 
which they perform spin measurement. The matrix elements has been obtained 
by a sample of 500 simulated data (From Ref. |231)- 

Summarizing, the ML technique can be used to estimate the density matrix 
of a quantum system. With respect to conventional quantum tomography this 
method has the great advantage of needing much smaller experimental samples, 
making experiments with low data rates feasible, however with a truncation of 
the Hilbert space dimension. We have shown that the method is general and the 
algorithm has solid methodological background, its reliability being confirmed 




(8.16) 
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in a number of Monte Carlo simulations. However, for increasing dimension of 
Hilbert spaces the method has exponential complexity. 



8.3 Gaussian-state estimation 

In this Section we review the ML determination method of Ref. jHS] for the 
parameters of Gaussian states. Such states represent the wide class of coherent, 
squeezed and thermal states, all of them being characterized by a Gaussian 
Wigner function. Apart from an irrelevant phase, we consider Wigncr functions 
of the form 

W{x, y) = — exp [e-2-(a; _ Re;,)2 + e'r^y _ Im^xf] ] , (8.17) 

and the ML technique with homodyne detection is applied to estimate the four 
real parameters A, r, Re/i and Im/x. The four parameters provide the number of 
thermal, squeezing and coherent-signal photons in the quantum state as follows 

Usq = sinh^ r , 

ncoh = Ia^I^ • (8.18) 
The density matrix p corresponding to the Wigner function in Eq. (|8.17|l writes 

p = D{^^) Sir) (^^) " " SHr) D\^,) , (8.19) 

nth + 1 \nth + 1/ 

where 5'(r) = exp[r(a2— a^2)/2] and _D(/i) = cx^{p,a) ~ p* a) denote the squeezing 
and displacement operators, respectively. 

The theoretical homodyne probability distribution at phase with respect 
to the local oscillator can be evaluated using Eq. (|2.7|) . and is given by the 
Gaussian 



2A2 



7r(e2'' cos^ + e -^"^ sit? (p) 



X exp^- ,^ ^ . , [x-Re{pe-^n] j • (8-20) 



e^'' cos^ If + e '^^ sin^ (p 



The log-likelihood function (|8.3|l for a set of N homodyne outcomes Xi at random 
phase ipi then writes as follows 



N 



2A2 



L ~ - log 2 

^ 2 7r(e2'' cos^ pi + e^^r gjj^ ^p.^ 



2A2 9 

[xi-Re{pe~"^^)] ■ (8.21) 



e^'' cos^ (pi + e 2'' sin^ (p. 
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The ML estimators A^; , fmi , Re/x^; and Imfimi are found upon maximizing Eq. 
H8.21|l versus A, r, Refi and Im^. 

In order to evaluate globally the state reconstruction, one considers the 
normalized overlap O between the theoretical and the estimated state 



'ml\ 



VTr[p2]Tr[p^ 



i.22) 



Notice that O = 1 iff p = p,„;. Through Monte-Carlo simulations, one always 
finds a value around unity, typically with statistical fluctuations over the third 
digit, for number of data samples N = 50000, quantum efficiency at homodyne 
detectors rj — 80%, and state parameters with the following ranges: nth < 3, 
ncoh < 5, and Usq < 3. Also with such a small number of data samples, the 
quality of the state reconstruction is so good that other physical quantities that 
are theoretically evaluated from the experimental values of Ami, fmh R^fJ-mi and 
Im/im/ are inferred very precisely. For example, in Ref. '33' the photon number 
probability of a squeezed thermal state has been evaluated, which is given by 
the integral 



{n\p\n) 



2tt C{(t),nth,rY 



+1 



.23) 



with C{(f>,nth,r) = {nth + ^)(e" 



2r , 



The comparison 



of the theoretical and the experimental results for a state with nth = 0.1 and 
nsq = 3 is reported in Fig. 18.41 The statistical error of the reconstructed number 
probability affects the third decimal digit, and is not visible on the scale of the 
plot. 

P(n) 



0.4 
0.3 
0.2 
0.1 



I 



.1.1.. 



01234567 



Figure 8.4: Photon-number probability of a squeezed-thermal state (thermal 
photons nth = 0.1, squeezing photons nsq = 3). Compare the reconstructed 
probabilities by means of the maximum likelihood method and homodyne de- 
tection (gray histogram) with the theoretical values (black histogram). Number 
of data samples N — 50000, quantum efficiency rj — 80%. The statistical error 
affects the third decimal digit, and it is not visible on the scale of the plot (From 
Ref. ES!)- 
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The estimation of parameters of Gaussian Wigner functions through the ML 
method allows one to estimate the parameters in quadratic Hamiltonians of the 
generic form 

H = aa + a*a^ + ipa^ a + ^£.a^ + ^Ca^^ ■ (8.24) 

In fact, the unitary evolution operator U — e~*^* preserves the Gaussian form 
of an input state with Gaussian Wigner function. In other words, one can use 
a known Gaussian state to probe and characterize an optical device described 
by a Hamiltonian as in Eq. H8.24|l . Assuming t = 1 without loss of generality, 
the Heisenberg evolution of the radiation mode a is given by 

C/^aC/ = 7a + 5a^ +/i , (8.25) 



with 



7 = cos( V'/'^ - m - *^==^ sH^^^^) , (8.26) 
S = -i—J= sin(V<^2 _ |^|2) 

ipa* — 



M = ^^^-|^(cos(y^^^) - 1) - *^^f=^ MV^P^) . 

For an input state p with known Wigner function Wp{p ,/?*), the corresponding 
output Wigner function writes 

Wupw (/3 , - Wp[i(3 ~ m)7* - (/3* " M*)'^ , {P* ~ ~ - f^)S*] ■ (8.27) 

Hence, by estimating the parameters j,S,fi and inverting Eqs. H8.26|l . one ob- 
tains the ML values for a, (p, and <^ of the Hamiltonian in Eq. (|8.24(l . The present 
example can be used in practical applications for the estimation of the gain of 
a phase-sensitive amplifier or equivalently to estimate a squeezing parameter. 
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Chapter 9 



Classical 
quantum 



imaging by 
tomography 



As we showed in Chapter 2, the development of quantum tomography has its 
origin in the inadequacy of classical imaging procedures to face the quantum 
problem of Wigner function reconstruction. In this chapter we briefly illustrate 
how to go back to classical imaging and profitably use quantum tomography as 
a tool for image reconstruction and compression: this is the method of Fictitious 
Photons Tomography of Ref . '34' . 

The problem of tomographic imaging is to recover a mass distribution m{x,y) 
in a 2-d slab from a finite collection of one dimensional projections. The situa- 
tion is schematically sketched in Fig. 19.1b where m{x, y) describes two circular 
holes in a uniform background. The tomographic machine, say an X-ray equip- 
ment, collects many stripe photos of the sample from various directions 6, and 
then numerically performs a mathematical transformation — the so called in- 
verse Radon transform jll6) — in order to reconstruct m(x, y) from its radial 
profiles at different 6''s. The problem which is of interest for us is when the 
radial profiles are not well-resolved digitalized functions, but actually represent 
the density distribution of random points, as if in our X-ray machine the beam 
is so weak that radial photos are just the collection of many small spots, each 
from a single X-ray photon (this situation is sketched in Fig. 19.1b ). It is obvious 
that this case can be reduced to the previous one by counting all points falling in 
a predetermined 1-d mesh, and giving radial profiles in form of histograms (this 
is what actually happens in a real machine, using arrays of photodetectors) . 
However, we want to use the whole available information from each "event" — 
i.e. the exact 1-d location of each spot — in a way which is independent on any 
predetermined mesh. In practice, this situation occurs when the signal is so 
weak and the machine resolution is so high (i.e. the mesh-step is so tiny) that 
only zero or one photon at most can be collected in each channel. As we will 
see, this low-signal/high-resolution case naturally brings the imaging problem 
into the domain of quantum tomography. Images are identified with Wigner 
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functions, such to obtain a description in terms of density matrices. These are 
stiU trace-class matrices (corresponding to " normahzable" images), but are no 
longer positive definite, because an "image" generally is not a genuine Wigner 
function and violates the Heisenberg relations on the complex plane (the phase 
space of a single mode of radiation). Hence, such density matrices are unphysi- 
cal: they are just a mathematical tool for imaging. This is the reason why this 
method has been named Fictitious Photons Tomography '34' . As we will see in 
the following, the image resolution improves by increasing the rank of the den- 
sity matrix, and in this way the present method also provides a new algorithm 
for image compression, which is suited to angular image scanning. 




Figure 9.1: (a) Tomography of a simple object: analytical transmission profiles 
are reported for 9 = 0, 7r/2. (b) The same case of (a), but for very weak signals: 
in this case the transmission profiles are given in terms of random points on a 
photographic plate (here obtained from a Monte Carlo simulation) (From Ref. 



9.1 From classical to quantum imaging 

We adopt the complex notation, with a = x + iy representing a point in the 
image plane. In this way a and a* are considered as independent variables, 
and the 2-d image — here denoted by the same symbol W{a,a*) used for the 
Wigner function — is just a generic real function of the point in the plane. In the 
most general situation W{a,a*) is defined on the whole complex plane, where 
it is normalized to some finite constant, and it is bounded from both below and 
above, with range representing the darkness nuance. For X-ray tomography 
W{a, a*) roughly represents the absorption coefficient as a function of the point 
a. We consider a linear absorption regime, i.e. the image extension is negligible 
with respect to the radiation absorption length in the medium. At the same 
time we neglect any diffraction effect. 
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As shown in Sec. 13.21 the customary imaging technique is based on the 
inverse Radon transform. A tomography of a two dimensional image W{a,a*) 
is a collection of one dimensional projections p{x, 9) at different values of the 
observation angle 9. We rewrite here the definition of the Radon transform of 
Wia,a*) 

p{x, 9) = r°^'^ W ((x + ^2/)e'^ (x - iy)e-^0) . (9.1) 

J-oo 

In Eq. 1)9. l|l x is the current coordinate along the direction orthogonal to the 
projection and y is the coordinate along the projection direction. The situation 
is depicted in Fig. 19. II where W{a,a.*) is plotted along with its p{x,9) profiles 
for 9 = 0,7r/2 for a couple of identical circular holes that are symmetrically 
disposed with respect to the origin. 

The reconstruction of the image VF(q;, a*) from its projections p{x, 9) — also 
called "back projection" — is given by the inverse Radon transform, which, fol- 
lowing the derivation in Sec. 13.21 leads to the filtering procedure 

r^Ppd^M^ii^, (9.2) 
Jo 27r X- ae 

where V denotes the Cauchy principal value and ag — Re(ae^'^). Eq. (|9.2|l 
is commonly used in conventional tomographic imaging (see, for example, Ref. 

inzi)- 

Let us now critically consider the above procedure in the case of very weak 
signals, namely when p(x, 9) just represents the probability distribution of ran- 
dom X-ray spots on a fine-mesh multichannel: this situation is sketched in Fig. 
19.1b . From Eq. H9.2|) one can recover W{a,a*) only when the analytical form 
of p{x, 9) is known. But the experimental outcomes of each projection actually 
are random data distributed according to p{x,9), whereas in order to recover 
W{a, a*) from Eq. I|9.2|l one has to evaluate the first order derivatives of p{x, 9). 
The need of the analytical form for projections p(x, 9) requires a filtering pro- 
cedure on data, usually obtained by "spHning" data in order to use Eq. H9.2|l . 

The above procedure leads to approximate image reconstructions, and the 
choice of any kind of smoothing parameter unavoidably affects in a systematic 
way the statistics of errors. In the following we show how quantum tomography 
can be used for conventional imaging in presence of weak signals, providing both 
ideally controlled resolution and reliable error statistics. 

The basic formula we will use is the expansion of the Wigner function in the 
number representation of Eqs. H6.22|l and H6.23|l . In practice, the Hilbert space 
has to be truncated at some finite dimension d^, and this sets the resolution 
for the reconstruction of W(a^a*). However, as we will show, this resolution 
can be chosen at will, independently on the number of experimental data. 

As previously noticed, in general an image does not correspond to a Wigner 
function of a physical state, due to the fact that the Heisenberg relations un- 
avoidably produce only smooth Wigner functions, whereas a conventional image 
can have very sharp edges. However, if one allows the density matrix to be no 
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longer positive definite (but still trace class), a correspondence with images is 
obtained, which holds in general. In this way every image is stored into a trace- 
class matrix pn,m via quantum tomography, and a convenient truncation of the 
matrix dimension d-u can be chosen. 



Symmetry 


p{x,9) 


p 


Isotropy 
X-axis mirror 
Y-axis mirror 
Inversion through the origin 


p{x,0) =p{x) 

p{x,Tr-e) =p{-x,e) 
p{x,Tr-e) ^p{x,e) 
p{x,e) =p{-x,0) 


Pn,m — Pn,n^n^ni 
Pn,m € K. 
iPn.m e R 
Pn,n+2d+l — 



Table 9.1: Geometrical symmetries of an image, analytical properties of pro- 
jections and algebraic properties of the corresponding matrix (From Ref. [SI)- 



The connection between images and matrices is the main point of this ap- 
proach: the information needed to reconstruct the image is stored in a d-^^ x d-u 
matrix. For suitably chosen dimension d-u the present method can also provide 
a procedure for image compression. Notice that the correspondence between 
images and trace-class matrices retains some symmetries of the image, which 
manifest as algebraic properties of the matrix Pn,m- For example an isotropic 
image (like a uniform circle centered at the origin) is stored in a diagonal matrix. 
Other symmetries are given in Tab. 19.11 

The truncated Hilbert space dimension dn sets the imaging resolution. The 
kind of resolution can be understood by studying the behavior of the kernels 
TZ[\n + d){n\]{x,9) ofEq. H3.44|l . which are averaged over the experimental data 
in order to obtain the matrix elements Pn,n+d- Outside a region that is almost 
independent of n and d, all functions TZ[\n + d) {n\]{x, 9) decrease exponentially, 
whereas inside this region they oscillate with a number of oscillations linearly 
increasing with 2n+d. This behavior produces the effects illustrated in Fig. 19.21 
where we report the tomographic reconstruction of the font "a" for increasing 
dimension dfi . The plot is obtained by numerically integrating the kernel func- 
tions from given analytic transmission profiles p{x, 9). As we see from Fig. 19.21 
both the radial and the angular resolutions improve versus d-u, making the de- 
tails of the image sharper and sharper already form a relatively small truncation 
d-H = 48. 

A quantitative measure of the precision of the tomographic reconstruction 
can be given in terms of the distance D between the true and the reconstructed 
image, which, in turn, coincides with the Hilbert distance D between the cor- 
responding density matrices. One has 

D = J d^a\AWia,a*)\^ :^Tr{Apy 

OO OQ OQ 

= 5^Ap^^„ + 2^^|Ap2„^,f , (9.3) 

n=Q n=OA=l 

where A[. . .] = [. . .]triie - [■ ■ ■]reconstructed- The Convergence of D versus dn 
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is given in Fig. 19.31 for a solid circle of unit radius centered at the origin. In 
this case the obtained density matrix has only diagonal elements, according to 
Tab. 19.11 These are given by 



where $(a,/3, z) denotes the confluent hypergeometric function of argument z 
and parameters a and /3. 

Insofar we have analyzed the method only on the basis of given analytic pro- 
files ■p{x, 9). As already said, however, the method is particularly advantageous 
in the weak-signal/high- resolution situation, where the imaging can be achieved 
directly from averaging the kernel functions on data. In this case the procedure 
allows to exploit the whole available experimental resolution, whereas the image 
resolution is set at will. In Fig. 19.41 we report a Monte Carlo simulation of an 
experimental tomographic reconstruction of the font "a" for increasing number 
of data. All plots are obtained at the maximum available dimension d-j-i = 48, 
and using F — 100 scanning phases. The situation occurring for small num- 
bers of data is given in the first plot, where the the highly resolved image still 
exhibits the natural statistical fluctuations due to the limited number of data. 
For larger sample the image appears sharper from the random background, and 
it is clearly recognizable for a number of data equal to 10^. The method is ef- 
ficient also from the computational point of view, as the time needed for image 
reconstruction is quadratic in the number of elements of the density matrix, and 
linear in the number of experimental data. Needless to say, imaging by quantum 
homodyne tomography is at the very early stages and further investigation is in 
order. 




(9.4) 
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f € 










Figure 9.2: Tomographic reconstruction of the font "a" for increasing dimension 
of the truncated matrix, dn ~ 2, 4, 8, 16, 32, 48. The plot is obtained by averag- 
ing the kernel function TZ[\n + o?)(n|](a::, 9) of Eq. (|3.44|) with assigned analytic 
transmission profiles p{x,9), and then using Eqs. H6.22() and (|6.22fl (From Ref. 

El). 
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Figure 9.3: Convergence of both trace and Hilbert distance D in Eq. (|9.3|l 
versus the dimensional truncation dn of the Hilbert space. Here the image is a 
uniform circle of unit radius centered at the origin. The reconstructed matrix 
elements are obtained as in Fig. 19.21 whereas the exact matrix elements are 
provided by Eq. (^31) (From Ref. 




Figure 9.4: Monte Carlo simulation of an experimental tomographic reconstruc- 
tion of the font "a" . The truncation dimension is fixed at dn = 48, and the 
number of scanning phases is F = 100. The plots correspond to 10^, 10^, 10^, 10^ 
data for each phase respectively (From Ref. |34|)- 
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